#include "basis_state.h"
//~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~start
// Jacobi~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~//
void Jacobi::build_JT(int n12max_t, int n3max_t, int l12max_t, int l3max_t) {
  n12max = n12max_t;
  n3max = n3max_t;
  l12max = l12max_t;
  l3max = l3max_t;
  Emax = 2 * (n12max + n3max) + l12max + l3max;

  J2max = (l12max + l3max) * 2 + 3;
  JT_num = (J2max + 1);
  state.resize(JT_num);
  int T3_2 = 1;
  // num=0;
  int JT_channel = -1;
  State_Jacobi state_t;
  for (int T_2 = 1; T_2 <= 3; T_2 = T_2 + 2) {
    for (int J_2 = 1; J_2 <= J2max; J_2 = J_2 + 2) {
      JT_channel++;
      // cout<<JT_channel<<endl;
      for (int T12_2 = 0; T12_2 <= 2; T12_2 = T12_2 + 2) {
        if (this->Tri(T_2, T12_2, T3_2) != 1)
          continue;
        for (int S12_2 = 0; S12_2 <= 2; S12_2 = S12_2 + 2) {
          for (int l12 = 0; l12 <= l12max; l12++) {
            int J12_2min = std::abs(l12 * 2 - S12_2);
            int J12_2max = l12 * 2 + S12_2;
            for (int J12_2 = J12_2min; J12_2 <= J12_2max; J12_2 = J12_2 + 2) {
              for (int l3 = 0; l3 <= l3max; l3++) {
                int j3_2min = std::abs(l3 * 2 - 1);
                int j3_2max = l3 * 2 + 1;
                for (int j3_2 = j3_2min; j3_2 <= j3_2max; j3_2 = j3_2 + 2) {
                  if (this->Tri(J12_2, j3_2, J_2) != 1)
                    continue;
                  int phase_t = std::abs((S12_2 + T12_2) / 2 + l12);
                  if ((phase_t % 2) != 1)
                    continue;
                  for (int n12 = 0; n12 <= n12max; n12++) {
                    for (int n3 = 0; n3 <= n3max; n3++) {
                      int E_x = 2 * (n12 + n3) + l12 + l3;
                      state_t.n12 = n12;
                      state_t.n3 = n3;
                      state_t.l12 = l12;
                      state_t.l3 = l3;
                      state_t.S12_2 = S12_2;
                      state_t.J12_2 = J12_2;
                      state_t.j3_2 = j3_2;
                      state_t.J_2 = J_2;
                      state_t.T_2 = T_2;
                      state_t.T12_2 = T12_2;
                      state_t.E = 2 * (n12 + n3) + l12 + l3;
                      state_t.Par = this->phase(l12 + l3);
                      state_t.index = -1;
                      state_t.flag = -1;
                      state_t.F = -1;
                      state_t.I = -1;
                      state_t.JT_channel = JT_channel;
                      //state_t.flag = T_2*40000+T12_2*8000+J_2*400+J12_2*20+j3_2 \
		      //+S12_2*0.1+l12*0.001+l3*0.00001+n12*0.000001+n3*0.0000001;
                      this->state[JT_channel].push_back(state_t);
                    }
                  }
                }
              }
            }
          }
        }
      }
    }
  }
  JT_channel++;
  channel.resize(JT_channel);
  for (int i = 0; i < JT_channel; i++) {
    if (state[i].size() == 0) {
      channel[i].num = 0;
      channel[i].J_2 = 0;
      channel[i].T_2 = 0;
      channel[i].E_min = 0;
      channel[i].E_max = 0;
    } else {
      channel[i].num = state[i].size();
      channel[i].J_2 = state[i][0].J_2;
      channel[i].T_2 = state[i][0].T_2;
      int E_min = state[i][0].E;
      int E_max = state[i][0].E;
      for (int j = 0; j < state[i].size(); j++) {
        if (E_min > state[i][j].E)
          E_min = state[i][j].E;
        if (E_max < state[i][j].E)
          E_max = state[i][j].E;
      }
      channel[i].E_min = E_min;
      channel[i].E_max = E_max;
    }
  }
  // calulate JT2channel
  JT2channel.resize(2); // 0: T=0.5; 1:T=1.5
  for (int i = 0; i <= 1; i++) {
    int T_2 = i * 2 + 1;
    JT2channel[i].resize(J2max + 1);
    for (int j = 0; j < JT2channel[i].size(); j++)
      JT2channel[i][j] = -1;
    for (int J2 = 1; J2 <= J2max; J2 = J2 + 2) {
      for (int k = 0; k < channel.size(); k++) {
        if (channel[k].T_2 == T_2 && channel[k].J_2 == J2) {
          JT2channel[i][J2] = k;
          break;
        }
      }
    }
  }
}

void Jacobi::build_JT(int Emax_t, int l12max_t, int l3max_t) {
  Emax = Emax_t;
  // n12max=std::ceil(Emax/2.0);
  // n3max=std::ceil(Emax/2.0);
  l12max = l12max_t;
  l3max = l3max_t;
  if (Emax - l12max - l3max < 0) {
    cout << "Error in the setting of Emax, Emax < l12max+l3max" << endl;
    exit(0);
  }

  J2max = (l12max + l3max) * 2 + 3;
  JT_num = (J2max + 1);
  state.resize(JT_num);
  int T3_2 = 1;
  // num=0;
  int JT_channel = -1;
  State_Jacobi state_t;
  for (int T_2 = 1; T_2 <= 3; T_2 = T_2 + 2) {
    for (int J_2 = 1; J_2 <= J2max; J_2 = J_2 + 2) {
      JT_channel++;
      cout << JT_channel << endl;
      for (int T12_2 = 0; T12_2 <= 2; T12_2 = T12_2 + 2) {
        if (this->Tri(T_2, T12_2, T3_2) != 1)
          continue;
        for (int S12_2 = 0; S12_2 <= 2; S12_2 = S12_2 + 2) {
          for (int l12 = 0; l12 <= l12max; l12++) {
            int J12_2min = std::abs(l12 * 2 - S12_2);
            int J12_2max = l12 * 2 + S12_2;
            for (int J12_2 = J12_2min; J12_2 <= J12_2max; J12_2 = J12_2 + 2) {
              for (int l3 = 0; l3 <= l3max; l3++) {
                int j3_2min = std::abs(l3 * 2 - 1);
                int j3_2max = l3 * 2 + 1;
                for (int j3_2 = j3_2min; j3_2 <= j3_2max; j3_2 = j3_2 + 2) {
                  if (this->Tri(J12_2, j3_2, J_2) != 1)
                    continue;
                  int phase_t = std::abs((S12_2 + T12_2) / 2 + l12);
                  if ((phase_t % 2) != 1)
                    continue;
                  int n12_max = std::ceil((Emax - l12 - l3) / 2.0);
                  for (int n12 = 0; n12 <= n12_max; n12++) {
                    int n3_max = std::ceil((Emax - l12 - l3 - 2 * n12) / 2.0);
                    for (int n3 = 0; n3 <= n3_max; n3++) {
                      int E_x = 2 * (n12 + n3) + l12 + l3;
                      if (E_x > Emax)
                        continue;
                      state_t.n12 = n12;
                      state_t.n3 = n3;
                      state_t.l12 = l12;
                      state_t.l3 = l3;
                      state_t.S12_2 = S12_2;
                      state_t.J12_2 = J12_2;
                      state_t.j3_2 = j3_2;
                      state_t.J_2 = J_2;
                      state_t.T_2 = T_2;
                      state_t.T12_2 = T12_2;
                      state_t.E = 2 * (n12 + n3) + l12 + l3;
                      state_t.Par = this->phase(l12 + l3);
                      state_t.index = -1;
                      state_t.flag = -1;
                      state_t.F = -1;
                      state_t.I = -1;
                      state_t.JT_channel = JT_channel;
                      //state_t.flag = T_2*40000+T12_2*8000+J_2*400+J12_2*20+j3_2 \
		      //+S12_2*0.1+l12*0.001+l3*0.00001+n12*0.000001+n3*0.0000001;
                      this->state[JT_channel].push_back(state_t);
                    }
                  }
                }
              }
            }
          }
        }
      }
    }
  }
  JT_channel++;

  channel.resize(JT_channel);
  for (int i = 0; i < JT_channel; i++) {
    if (state[i].size() == 0) {
      channel[i].num = 0;
      channel[i].J_2 = 0;
      channel[i].T_2 = 0;
      channel[i].E_min = 0;
      channel[i].E_max = 0;
    } else {
      channel[i].num = state[i].size();
      channel[i].J_2 = state[i][0].J_2;
      channel[i].T_2 = state[i][0].T_2;
      int E_min = state[i][0].E;
      int E_max = state[i][0].E;
      for (int j = 0; j < state[i].size(); j++) {
        if (E_min > state[i][j].E)
          E_min = state[i][j].E;
        if (E_max < state[i][j].E)
          E_max = state[i][j].E;
      }
      channel[i].E_min = E_min;
      channel[i].E_max = E_max;
    }
  }
  // calulate JT2channel
  JT2channel.resize(2); // 0: T=0.5; 1:T=1.5
  for (int i = 0; i <= 1; i++) {
    int T_2 = i * 2 + 1;
    JT2channel[i].resize(J2max + 1);
    for (int j = 0; j < JT2channel[i].size(); j++)
      JT2channel[i][j] = -1;
    for (int J2 = 1; J2 <= J2max; J2 = J2 + 2) {
      for (int k = 0; k < channel.size(); k++) {
        if (channel[k].T_2 == T_2 && channel[k].J_2 == J2) {
          JT2channel[i][J2] = k;
          break;
        }
      }
    }
  } /**/
}

void Jacobi::build_JT(int Emax_t, int J2max_t) {
  Emax = Emax_t;
  J2max = J2max_t;
  // cout << "Jacobi::build_JT Emax =" << Emax << endl;
  // cout << "Jacobi::build_JT J2max =" << J2max << endl;
  // n12max=std::ceil(Emax/2.0);
  // n3max=std::ceil(Emax/2.0);
  // l12max=l12max_t;
  // l3max=l3max_t;
  // J2max = J2max_t + 3;

  // J2max=(Emax)*2+3;
  // J2max=(l12max+l3max)*2+3;
  JT_num = (J2max + 1);
  state.resize(JT_num);
  int T3_2 = 1;
  // num=0;
  int JT_channel = 0;
  State_Jacobi state_t;
  for (int T_2 = 1; T_2 <= 3; T_2 = T_2 + 2) {
    for (int J_2 = 1; J_2 <= J2max; J_2 = J_2 + 2) {
      // JT_channel++;
      // cout<<JT_channel<<endl;
      for (int T12_2 = 0; T12_2 <= 2; T12_2 = T12_2 + 2) {
        if (this->Tri(T_2, T12_2, T3_2) != 1)
          continue;
        for (int S12_2 = 0; S12_2 <= 2; S12_2 = S12_2 + 2) {
          int l12_max = Emax * 2;
          for (int l12 = 0; l12 <= l12_max; l12++) {
            int J12_2min = std::abs(l12 * 2 - S12_2);
            int J12_2max = l12 * 2 + S12_2;
            for (int J12_2 = J12_2min; J12_2 <= J12_2max; J12_2 = J12_2 + 2) {
              int l3_max = Emax * 2;
              for (int l3 = 0; l3 <= l3_max; l3++) {
                int j3_2min = std::abs(l3 * 2 - 1);
                int j3_2max = l3 * 2 + 1;
                for (int j3_2 = j3_2min; j3_2 <= j3_2max; j3_2 = j3_2 + 2) {
                  if (this->Tri(J12_2, j3_2, J_2) != 1)
                    continue;
                  int phase_t = std::abs((S12_2 + T12_2) / 2 + l12);
                  if ((phase_t % 2) != 1)
                    continue;
                  // int n12_max=std::ceil((Emax-l12-l3)/2.0)+1;
                  int n12_max = Emax;
                  for (int n12 = 0; n12 <= n12_max; n12++) {
                    // int n3_max=std::ceil((Emax-l12-l3-2*n12)/2.0)+1;
                    int n3_max = Emax;
                    for (int n3 = 0; n3 <= n3_max; n3++) {
                      int E_x = 2 * (n12 + n3) + l12 + l3;
                      if (E_x > Emax)
                        continue;
                      state_t.n12 = n12;
                      state_t.n3 = n3;
                      state_t.l12 = l12;
                      state_t.l3 = l3;
                      state_t.S12_2 = S12_2;
                      state_t.J12_2 = J12_2;
                      state_t.j3_2 = j3_2;
                      state_t.J_2 = J_2;
                      state_t.T_2 = T_2;
                      state_t.T12_2 = T12_2;
                      state_t.E = 2 * (n12 + n3) + l12 + l3;
                      state_t.Par = this->phase(l12 + l3);
                      state_t.index = -1;
                      state_t.flag = -1;
                      state_t.F = -1;
                      state_t.I = -1;
                      state_t.JT_channel = JT_channel;
                      //state_t.flag = T_2*40000+T12_2*8000+J_2*400+J12_2*20+j3_2 \
		      //+S12_2*0.1+l12*0.001+l3*0.00001+n12*0.000001+n3*0.0000001;
                      this->state[JT_channel].push_back(state_t);
                    }
                  }
                }
              }
            }
          }
        }
      }
      if (state[JT_channel].size() > 0)
        JT_channel++;
    }
  }
  // JT_channel++;

  channel.resize(JT_channel);
  for (int i = 0; i < JT_channel; i++) {
    if (state[i].size() == 0) {
      channel[i].num = 0;
      channel[i].J_2 = 0;
      channel[i].T_2 = 0;
      channel[i].E_min = 0;
      channel[i].E_max = 0;
    } else {
      channel[i].num = state[i].size();
      channel[i].J_2 = state[i][0].J_2;
      channel[i].T_2 = state[i][0].T_2;
      int E_min = state[i][0].E;
      int E_max = state[i][0].E;
      for (int j = 0; j < state[i].size(); j++) {
        if (E_min > state[i][j].E)
          E_min = state[i][j].E;
        if (E_max < state[i][j].E)
          E_max = state[i][j].E;
      }
      channel[i].E_min = E_min;
      channel[i].E_max = E_max;
    }
  }
  // calulate JT2channel
  JT2channel.resize(2); // 0: T=0.5; 1:T=1.5
  for (int i = 0; i <= 1; i++) {
    int T_2 = i * 2 + 1;
    JT2channel[i].resize(J2max + 1);
    for (int j = 0; j < JT2channel[i].size(); j++)
      JT2channel[i][j] = -1;
    for (int J2 = 1; J2 <= J2max; J2 = J2 + 2) {
      for (int k = 0; k < channel.size(); k++) {
        if (channel[k].T_2 == T_2 && channel[k].J_2 == J2) {
          JT2channel[i][J2] = k;
          break;
        }
      }
    }
  } /**/
}

void Jacobi::build_JT(int Emax_t) {
  Emax = Emax_t;
  // n12max=std::ceil(Emax/2.0);
  // n3max=std::ceil(Emax/2.0);
  // l12max=l12max_t;
  // l3max=l3max_t;

  J2max = (Emax)*2 + 3;
  // J2max=(l12max+l3max)*2+3;
  JT_num = (J2max + 1);
  state.resize(JT_num);
  int T3_2 = 1;
  // num=0;
  int JT_channel = 0;
  State_Jacobi state_t;
  for (int T_2 = 1; T_2 <= 3; T_2 = T_2 + 2) {
    for (int J_2 = 1; J_2 <= J2max; J_2 = J_2 + 2) {
      // JT_channel++;
      // cout<<JT_channel<<endl;
      for (int T12_2 = 0; T12_2 <= 2; T12_2 = T12_2 + 2) {
        if (this->Tri(T_2, T12_2, T3_2) != 1)
          continue;
        for (int S12_2 = 0; S12_2 <= 2; S12_2 = S12_2 + 2) {
          int l12_max = Emax;
          for (int l12 = 0; l12 <= l12_max; l12++) {
            int J12_2min = std::abs(l12 * 2 - S12_2);
            int J12_2max = l12 * 2 + S12_2;
            for (int J12_2 = J12_2min; J12_2 <= J12_2max; J12_2 = J12_2 + 2) {
              int l3_max = Emax - l12;
              for (int l3 = 0; l3 <= l3_max; l3++) {
                int j3_2min = std::abs(l3 * 2 - 1);
                int j3_2max = l3 * 2 + 1;
                for (int j3_2 = j3_2min; j3_2 <= j3_2max; j3_2 = j3_2 + 2) {
                  if (this->Tri(J12_2, j3_2, J_2) != 1)
                    continue;
                  int phase_t = std::abs((S12_2 + T12_2) / 2 + l12);
                  if ((phase_t % 2) != 1)
                    continue;
                  int n12_max = std::ceil((Emax - l12 - l3) / 2.0) + 1;
                  for (int n12 = 0; n12 <= n12_max; n12++) {
                    int n3_max =
                        std::ceil((Emax - l12 - l3 - 2 * n12) / 2.0) + 1;
                    for (int n3 = 0; n3 <= n3_max; n3++) {
                      int E_x = 2 * (n12 + n3) + l12 + l3;
                      if (E_x > Emax)
                        continue;
                      state_t.n12 = n12;
                      state_t.n3 = n3;
                      state_t.l12 = l12;
                      state_t.l3 = l3;
                      state_t.S12_2 = S12_2;
                      state_t.J12_2 = J12_2;
                      state_t.j3_2 = j3_2;
                      state_t.J_2 = J_2;
                      state_t.T_2 = T_2;
                      state_t.T12_2 = T12_2;
                      state_t.E = 2 * (n12 + n3) + l12 + l3;
                      state_t.Par = this->phase(l12 + l3);
                      state_t.index = -1;
                      state_t.flag = -1;
                      state_t.F = -1;
                      state_t.I = -1;
                      state_t.JT_channel = JT_channel;
                      //state_t.flag = T_2*40000+T12_2*8000+J_2*400+J12_2*20+j3_2 \
		      //+S12_2*0.1+l12*0.001+l3*0.00001+n12*0.000001+n3*0.0000001;
                      this->state[JT_channel].push_back(state_t);
                    }
                  }
                }
              }
            }
          }
        }
      }
      if (state[JT_channel].size() > 0)
        JT_channel++;
    }
  }
  // JT_channel++;
  channel.resize(JT_channel);
  for (int i = 0; i < JT_channel; i++) {
    if (state[i].size() == 0) {
      channel[i].num = 0;
      channel[i].J_2 = 0;
      channel[i].T_2 = 0;
      channel[i].E_min = 0;
      channel[i].E_max = 0;
    } else {
      channel[i].num = state[i].size();
      channel[i].J_2 = state[i][0].J_2;
      channel[i].T_2 = state[i][0].T_2;
      int E_min = state[i][0].E;
      int E_max = state[i][0].E;
      for (int j = 0; j < state[i].size(); j++) {
        if (E_min > state[i][j].E)
          E_min = state[i][j].E;
        if (E_max < state[i][j].E)
          E_max = state[i][j].E;
      }
      channel[i].E_min = E_min;
      channel[i].E_max = E_max;
    }
  }
  // calulate JT2channel
  JT2channel.resize(2); // 0: T=0.5; 1:T=1.5
  for (int i = 0; i <= 1; i++) {
    int T_2 = i * 2 + 1;
    JT2channel[i].resize(J2max + 1);
    for (int j = 0; j < JT2channel[i].size(); j++)
      JT2channel[i][j] = -1;
    for (int J2 = 1; J2 <= J2max; J2 = J2 + 2) {
      for (int k = 0; k < channel.size(); k++) {
        if (channel[k].T_2 == T_2 && channel[k].J_2 == J2) {
          JT2channel[i][J2] = k;
          break;
        }
      }
    }
  }
}

void Jacobi::print() {
  // cout<<"total number of basis_JT is : "<<endl;
  // cout<<num<<endl;
  cout << "number of JT channel is : " << endl;
  cout << JT_num << endl;
  for (int i = 0; i < JT_num; i++) {
    cout << "===JT channel : " << i << " ====" << endl;
    cout << "num: \t"
         << "2*T \t"
         << "2*T12 \t"
         << "2*J \t"
         << "2*J12 \t"
         << "2*j3 \t"
         << "2*S12 \t"
         << "l12 \t"
         << "l3 \t"
         << "n12 \t"
         << "n3 \t"
         << "E \t"
         << "index \t"
         << "JT_channel \t"
         << "flag" << endl;
    for (int j = 0; j < state[i].size(); j++) {
      cout << j + 1 << "\t";
      cout << state[i][j].T_2 << "\t" << state[i][j].T12_2 << "\t"
           << state[i][j].J_2 << "\t" << state[i][j].J12_2 << "\t";
      cout << state[i][j].j3_2 << "\t";
      cout << state[i][j].S12_2 << "\t" << state[i][j].l12 << "\t"
           << state[i][j].l3 << "\t" << state[i][j].n12 << "\t";
      cout << state[i][j].n3 << "\t" << state[i][j].E << "\t"
           << state[i][j].index << "\t";
      cout << state[i][j].JT_channel << "\t\t" << state[i][j].flag << endl;
    }
  }
}

void Jacobi::print(int channel_num) {
  // cout<<"total number of basis_JT is : "<<endl;
  // cout<<num<<endl;
  cout << "number of JT channel is : " << endl;
  cout << JT_num << endl;
  int i = channel_num;
  // for(int i=0;i<JT_num;i++){
  cout << "===JT channel : " << i << " ====" << endl;
  cout << "num: \t"
       << "2*T \t"
       << "2*T12 \t"
       << "2*J \t"
       << "2*J12 \t"
       << "2*j3 \t"
       << "2*S12 \t"
       << "l12 \t"
       << "l3 \t"
       << "n12 \t"
       << "n3 \t"
       << "E \t"
       << "index \t"
       << "JT_channel \t"
       << "flag" << endl;
  for (int j = 0; j < state[i].size(); j++) {
    cout << j + 1 << "\t";
    cout << state[i][j].T_2 << "\t" << state[i][j].T12_2 << "\t"
         << state[i][j].J_2 << "\t" << state[i][j].J12_2 << "\t";
    cout << state[i][j].j3_2 << "\t";
    cout << state[i][j].S12_2 << "\t" << state[i][j].l12 << "\t"
         << state[i][j].l3 << "\t" << state[i][j].n12 << "\t";
    cout << state[i][j].n3 << "\t" << state[i][j].E << "\t" << state[i][j].index
         << "\t";
    cout << state[i][j].JT_channel << "\t\t" << state[i][j].flag << endl;
  }
  //}
}

void Jacobi::print(std::string filename) {
  std::ofstream file;
  file.open(filename);
  //file<<"Jacobi basis with l3max = "<<l3max<<", \t l12max = "<<l12max<<", \t n12max = "<<n12max \
      //    <<"\t n3max = "<<n3max<<endl;
  file << "Energy cut-off is : " << Emax << " ( hbar-omega)" << endl;
  file << "total number of basis is : " << endl;
  int total_num = 0;
  for (int i = 0; i < state.size(); i++) {
    total_num += state[i].size();
  }
  file << total_num << endl;
  file << state.size() << "\t JT channel num" << endl;
  for (int i = 0; i < state.size(); i++) {
    file << state[i].size() << "  ";
  }
  file << endl;
  file << "num: \t"
       << "2*T \t"
       << "2*T12 \t"
       << "2*J \t"
       << "2*J12 \t"
       << "2*j3 \t"
       << "2*S12 \t"
       << "l12 \t"
       << "l3 \t"
       << "n12 \t"
       << "n3 \t"
       << "E \t"
       << "index\t "
       << "JT_channel\t F\t I" << endl;
  for (int i = 0; i < state.size(); i++) {
    file << "channel No. : " << i << endl;
    for (int j = 0; j < state[i].size(); j++) {
      file << j + 1 << "\t";
      file << state[i][j].T_2 << "\t" << state[i][j].T12_2 << "\t"
           << state[i][j].J_2 << "\t" << state[i][j].J12_2 << "\t";
      file << state[i][j].j3_2 << "\t";
      file << state[i][j].S12_2 << "\t" << state[i][j].l12 << "\t"
           << state[i][j].l3 << "\t" << state[i][j].n12 << "\t";
      file << state[i][j].n3 << "\t" << state[i][j].E << "\t"
           << state[i][j].index << "\t" << state[i][j].JT_channel
           << "\t "; // << "\t";
      // file << state[i][j].F << "\t";
      // file << state[i][j].I << "\t";

      if (state[i][j].F == 1)
        file << "*";
      file << "\t";
      if (state[i][j].I == 1)
        file << "* \t"; /**/
      file << endl;
    }
  }
  file.close();
}

void Jacobi::print_channel() {
  int size = this->channel.size();
  cout << " Jacobi channel Total Num: " << size << endl;
  cout << "index \t"
       << "num \t"
       << "J2 \t"
       << "T_2 \t"
       << "E_min \t"
       << "E_max \t" << endl;
  for (int i = 0; i < size; i++) {
    cout << i << "\t" << channel[i].num << "\t" << channel[i].J_2 << "\t"
         << channel[i].T_2 << "\t" << channel[i].E_min << "\t"
         << channel[i].E_max << endl;
  }
}
void Jacobi::print_channel(std::string filename) {
  std::ofstream file;
  file.open(filename);
  int size = this->channel.size();
  file << " Jacobi channel Total Num: " << size << endl;
  file << "index \t"
       << "num \t"
       << "J2 \t"
       << "T_2 \t"
       << "E_min \t"
       << "E_max \t" << endl;
  for (int i = 0; i < size; i++) {
    file << i << "\t" << channel[i].num << "\t" << channel[i].J_2 << "\t"
         << channel[i].T_2 << "\t" << channel[i].E_min << "\t"
         << channel[i].E_max << endl;
  }
  file << "E_min and E_max index:" << endl;
  file << "\t \t   "
       << "E : 0\t1\t2 \t3 \t4 \t5 ....." << endl;
  for (int i = 0; i < size; i++) {
    file << "channel : " << i << endl;
    file << "\t E_min index : ";
    for (int j = 0; j < channel[i].indexE_beg.size(); j++) {
      file << channel[i].indexE_beg[j] << "\t";
    }
    file << endl;
    file << "\t E_max index : ";
    for (int j = 0; j < channel[i].indexE.size(); j++) {
      file << channel[i].indexE[j] << "\t";
    }
    file << endl;
  }
  file.close();
}

void Jacobi::print_Jacconf(std::string filename) {
  std::ofstream file;
  // file.open("data/Jacconf.dat");
  int num = 0;
  for (int ch = 0; ch < Jacconf.size(); ch++) {
    num += Jacconf[ch].size();
  }
  file.open(filename);
  file << "Jacobi conf total : " << num << endl;
  for (int ch = 0; ch < Jacconf.size(); ch++) {
    file << "\t \t ch : " << ch
         << " Jacconf[ch].size() : " << Jacconf[ch].size() << endl;
  }
  for (int ch = 0; ch < Jacconf.size(); ch++) {
    file << "\t == channel : " << ch << "\t == " << endl;
    for (int i = 0; i < Jacconf[ch].size(); i++) {
      file << i << "\t" << Jacconf[ch][i].F << "\t" << Jacconf[ch][i].I << endl;
    }
  }
  file.close();
}

int Jacobi::Tri(int &j2a, int &j2b, int &j2c) {
  // I+J>=K, I+K>=J, J+K>=I,
  // I/2+J/2+K/2 = INTEGER.
  // TRI=1, WHEN TRIADIC CONDITION IS FULFILLED, TRI=-1 OTHERWISE
  int L2 = j2a + j2b + j2c;
  if (int(L2) != (L2 / 2) * 2)
    return -1;
  L2 = L2 / 2;
  if (j2a * j2b * j2c < 0)
    return -1;
  if ((L2 - j2a) * (L2 - j2b) * (L2 - j2c) < 0)
    return -1;
  return 1;
}

int Jacobi::Tri_state(State_Jacobi &state) {
  int unit = 1;
  int l12_2 = 2 * state.l12;
  int l3_2 = 2 * state.l3;
  if (Tri(l12_2, state.S12_2, state.J12_2) != 1) {
    return -1;
  }
  if (Tri(l3_2, unit, state.j3_2) != 1) {
    return -2;
  }
  if (Tri(state.J12_2, state.j3_2, state.J_2) != 1) {
    return -3;
  }
  if (Tri(state.T12_2, unit, state.T_2) != 1) {
    return -4;
  }
  return 1;
}

void Jacobi::check() {
  for (int i = 0; i < state.size(); i++) {
    for (int j = 0; j < state[i].size(); j++) {
      if (Tri_state(state[i][j]) < 0) {
        cout << " !!! wrong basis :" << i
             << "  wrong type : " << Tri_state(state[i][j]) << endl;
        cout << state[i][j].l3 << "\t" << state[i][j].j3_2 << endl;
      }
    }
  }
}

void Jacobi::select() {

  // select good by flag == 1
  // cout << "Start to select Jacobi state" << endl;
  for (int i = 0; i < state.size(); i++) {
    auto it = state[i].begin();
    while (it != state[i].end()) {
      if (it->flag == -1) {
        it = state[i].erase(it);
      } else {
        it++;
      }
    }
  }
  this->sort();
  // cout << "End to selete Jacobi state" << endl;

  /*// change the channel
  for (int i = 0; i < channel.size(); i++) {
    if (state[i].size() == 0) {
      channel[i].num = 0;
      // channel[i].J_2 = 0;
      // channel[i].T_2 = 0;
      channel[i].E_min = 0;
      channel[i].E_max = 0;
    } else {
      channel[i].num = state[i].size();
      channel[i].J_2 = state[i][0].J_2;
      channel[i].T_2 = state[i][0].T_2;
      int E_min = state[i][0].E;
      int E_max = state[i][0].E;
      for (int j = 0; j < state[i].size(); j++) {
        if (E_min >= state[i][j].E)
          E_min = state[i][j].E;
        if (E_max <= state[i][j].E)
          E_max = state[i][j].E;
      }
      channel[i].E_min = E_min;
      channel[i].E_max = E_max;
    }
  }
  // mark the index
  for (int i = 0; i < channel.size(); i++) {
    int E_min = channel[i].E_min;
    int E_max = channel[i].E_max;
    channel[i].indexE.resize(E_max + 1);
    channel[i].indexE_beg.resize(E_max + 1);
    for (auto it = channel[i].indexE.begin(); it != channel[i].indexE.end();
         it++)
      channel[i].indexE[*it] = 0;
    for (auto it = channel[i].indexE_beg.begin();
         it != channel[i].indexE_beg.end(); it++)
      channel[i].indexE_beg[*it] = 0;
    int index = 0;
    for (int k = 0; k < state[i].size(); k++) {
      state[i][k].index = index;
      int E = state[i][k].E;
      channel[i].indexE[E] = index;
      index++;
    }
    // indexE_beg
    for (int E = E_min; E <= E_max; E++) {
      int indexE_beg = 0;
      for (int k = 0; k < state[i].size(); k++) {
        if (state[i][k].E == E) {
          indexE_beg = k;
          break;
        }
      }
      channel[i].indexE_beg[E] = indexE_beg;
    }
  }/**/
}

void Jacobi::select_FI_T(vector<vector<State_TPconf_E>> TP_conf_E_tot) {
  int E_F_max = TP_conf_E_tot[0][0].E_F;
  int E_I_max = TP_conf_E_tot[0][0].E_I;
  for (int k = 0; k < TP_conf_E_tot.size(); k++) {
    for (int i = 0; i < TP_conf_E_tot[k].size(); i++) {
      if (E_F_max < TP_conf_E_tot[k][i].E_F)
        E_F_max = TP_conf_E_tot[k][i].E_F;
      if (E_I_max < TP_conf_E_tot[k][i].E_I)
        E_I_max = TP_conf_E_tot[k][i].E_I;
    }

    for (int ch = 0; ch < state.size(); ch++) {
      for (int i = 0; i < state[ch].size(); i++) {
        int E = state[ch][i].E;
        if (E <= E_F_max)
          state[ch][i].F = 1;
        if (E <= E_I_max)
          state[ch][i].I = 1;
      }
    }
  }
  // !!! Delete the state
  /*
  for (int ch = 0; ch < state.size(); ch++) {
    auto it = state[ch].begin();
    while (it != state[ch].end()) {
      if (it->F < 0 && it->I < 0) {
        it = state[ch].erase(it);
      } else {
        it++;
      }
    }
  }/**/
}
void Jacobi::Delet_Sort_T(TripleP &TP) {
  // call Select_FI_T first !!!
  // 1. Delete by .F = -1 && .I = -1
  //    Delete E > E_TP for each channel
  // 2. Sort for index and channel

  // check Emax of each channel
  /*for (int i = 0; i < channel.size(); i++) {
    int J2 = channel[i].J_2;
    int T2 = channel[i].T_2;
    cout << " ## J2 = " << J2 << "\t T2 = " << T2 << endl;
    int ch_TP = TP.FindC(J2, T2);
    cout << "TP channel : " << ch_TP << endl;
    int E_TP = TP.channel[ch_TP].E_max;

    if (channel[i].E_max < E_TP) {
      cout << " Jacobi Emax < Emax of TP in channel : " << i << endl;
      cout << " Please change a bigger Emax input in Jacobi.build_JT(Emax, "
              "J2max)"
           << endl;
      cout << "channel[i].E_max = " << channel[i].E_max << "\t E_TP = " << E_TP
           << endl;
      exit(0);
    }
  }/**/
  // delete state which E > E_TP
  vector<vector<State_Jacobi>> state_t;
  int JT_channel = state.size();
  state_t = state;
  state.clear();
  state.resize(JT_channel);

  int Emax_t = -1;
  for (int ch = 0; ch < JT_channel; ch++) {
    for (int i = 0; i < state_t[ch].size(); i++) {
      if (state_t[ch][i].F < 0 && state_t[ch][i].I < 0)
        continue;

      state[ch].push_back(state_t[ch][i]);
      if (state_t[ch][i].E > Emax_t)
        Emax_t = state_t[ch][i].E;
    }
  }
  this->Emax = Emax_t;

  auto it = state.begin();
  while (it != state.end()) {
    if (it->size() == 0) {
      it = state.erase(it);
    } else {
      it++;
    }
  }
}

void Jacobi::build_conf_T(vector<vector<State_TPconf_E>> TP_conf_E_tot) {

  int size_ch = this->state.size();
  this->Jacconf.resize(size_ch);
  int num = 0;
  for (int ch = 0; ch < size_ch; ch++) {
    for (int f = 0; f < state[ch].size(); f++) {
      for (int i = 0; i < state[ch].size(); i++) {
        int flag_conf = -1;
        for (int k = 0; k < TP_conf_E_tot.size(); k++) {
          for (int j = 0; j < TP_conf_E_tot[k].size(); j++) {
            int E_F = TP_conf_E_tot[k][j].E_F;
            int E_I = TP_conf_E_tot[k][j].E_I;
            int flag = -1;
            int E_diff = 0;
            if (E_F <= E_I) {
              flag = 1;
              E_diff = E_I - E_F;
            } else {
              flag = 2;
              E_diff = E_F - E_I;
            }
            int E_min_max = std::min(E_F, E_I);
            for (int E_min = 0; E_min <= E_min_max; E_min++) {
              int E_F_jac = -1;
              int E_I_jac = -1;
              if (flag == 1) {
                E_F_jac = E_min;
                E_I_jac = E_min + E_diff;
              } else {
                E_F_jac = E_min + E_diff;
                E_I_jac = E_min;
              }

              if (state[ch][f].E == E_F_jac && state[ch][i].E == E_I_jac) {
                flag_conf = 1;
                break;
              }
            }
            if (flag_conf == 1)
              break;
          }
          if (flag_conf == 1)
            break;
        }
        if (flag_conf == 1) {
          State_Conf conf_t;
          conf_t.F = f;
          conf_t.I = i;
          Jacconf[ch].push_back(conf_t);
          num++;
        }
      }
    }
  }
  /*std::ofstream file;
  file.open("data/Jacconf.dat");
  file << "Jacobi conf total : " << num << endl;
  for (int ch = 0; ch < Jacconf.size(); ch++) {
    file << "\t == channel : " << ch << "\t == " << endl;
    for (int i = 0; i < Jacconf[ch].size(); i++) {
      file << i << "\t" << Jacconf[ch][i].F << "\t" << Jacconf[ch][i].I << endl;
    }
  }
  file.close();/**/
}
void Jacobi::sort() {
  // sort by E
  auto it = state.begin();
  while (it != state.end()) {
    if (it->size() == 0) {
      it = state.erase(it);
    } else {
      it++;
    }
  }
  int chan_size = state.size();
  for (int i = 0; i < chan_size; i++)
    stable_sort(this->state[i].begin(), this->state[i].end(),
                [](State_Jacobi a, State_Jacobi b) {
                  return a.E * 100000 + a.T12_2 * 10000 + a.J12_2 * 100 +
                             a.j3_2 + 0.1 * a.S12_2 + 0.001 * a.l12 +
                             0.0001 * a.l3 + 0.00001 * a.n12 +
                             0.000001 * a.n3 <=
                         b.E * 100000 + b.T12_2 * 10000 + b.J12_2 * 100 +
                             b.j3_2 + 0.1 * b.S12_2 + 0.001 * b.l12 +
                             0.0001 * b.l3 + 0.00001 * b.n12 + 0.000001 * b.n3;
                });
  for (int ch = 0; ch < chan_size; ch++) {
    for (int i = 0; i < state[ch].size(); i++) {
      state[ch][i].index = i;
      state[ch][i].JT_channel = ch;
    }
  }
  channel.clear();
  channel.resize(chan_size);
  for (int i = 0; i < chan_size; i++) {
    if (state[i].size() == 0) {
      channel[i].num = 0;
      channel[i].J_2 = 0;
      channel[i].T_2 = 0;
      channel[i].E_min = 0;
      channel[i].E_max = 0;
    } else {
      channel[i].num = state[i].size();
      channel[i].J_2 = state[i][0].J_2;
      channel[i].T_2 = state[i][0].T_2;
      int E_min = state[i][0].E;
      int E_max = state[i][0].E;
      for (int j = 0; j < state[i].size(); j++) {
        if (E_min > state[i][j].E)
          E_min = state[i][j].E;
        if (E_max < state[i][j].E)
          E_max = state[i][j].E;
      }
      channel[i].E_min = E_min;
      channel[i].E_max = E_max;
    }
  }
  // calulate JT2channel
  JT2channel.resize(2); // 0: T=0.5; 1:T=1.5
  for (int i = 0; i <= 1; i++) {
    int T_2 = i * 2 + 1;
    JT2channel[i].resize(J2max + 1);
    for (int j = 0; j < JT2channel[i].size(); j++)
      JT2channel[i][j] = -1;
    for (int J2 = 1; J2 <= J2max; J2 = J2 + 2) {
      for (int k = 0; k < channel.size(); k++) {
        if (channel[k].T_2 == T_2 && channel[k].J_2 == J2) {
          JT2channel[i][J2] = k;
          break;
        }
      }
    }
  } /**/

  // mark the index
  for (int i = 0; i < chan_size; i++) {
    int E_min = channel[i].E_min;
    int E_max = channel[i].E_max;
    channel[i].indexE.resize(E_max + 1);
    channel[i].indexE_beg.resize(E_max + 1);
    for (auto it = channel[i].indexE.begin(); it != channel[i].indexE.end();
         it++)
      channel[i].indexE[*it] = 0;
    for (auto it = channel[i].indexE_beg.begin();
         it != channel[i].indexE_beg.end(); it++)
      channel[i].indexE_beg[*it] = 0;
    int index = 0;
    for (int k = 0; k < state[i].size(); k++) {
      state[i][k].index = index;
      int E = state[i][k].E;
      channel[i].indexE[E] = index;
      index++;
    }
    // indexE_beg
    for (int E = E_min; E <= E_max; E++) {
      int indexE_beg = 0;
      for (int k = 0; k < state[i].size(); k++) {
        if (state[i][k].E == E) {
          indexE_beg = k;
          break;
        }
      }
      channel[i].indexE_beg[E] = indexE_beg;
    }
  }
  // find l12max, l3max, n12max, n3max
  int l12max_t = 0;
  int l3max_t = 0;
  int l12min_t = 100;
  int l3min_t = 100;

  int n12max_t = 0;
  int n3max_t = 0;
  int n12min_t = 0;
  int n3min_t = 0;

  int J12_2max_t = 0;
  int j3_2max_t = 1;
  int all_size_t = 0;

  int E_max_t = 0;
  int E_min_t = 100;

  int l12max_F_t = 0;
  int l3max_F_t = 0;
  int l12min_F_t = 100;
  int l3min_F_t = 100;

  int l12max_I_t = 0;
  int l3max_I_t = 0;
  int l12min_I_t = 100;
  int l3min_I_t = 100;

  int E_F_max_t = 0;
  int E_I_max_t = 0;
  int E_F_min_t = 100;
  int E_I_min_t = 100;

  for (int ch = 0; ch < state.size(); ch++) {
    int size_i = state[ch].size();
    for (int i = 0; i < size_i; i++) {
      all_size_t++;

      if (state[ch][i].n12 > n12max_t)
        n12max_t = state[ch][i].n12;
      if (state[ch][i].n3 > n3max_t)
        n3max_t = state[ch][i].n3;

      if (state[ch][i].n12 < n12min_t)
        n12min_t = state[ch][i].n12;
      if (state[ch][i].n3 < n3min_t)
        n3min_t = state[ch][i].n3;

      if (state[ch][i].l12 > l12max_t)
        l12max_t = state[ch][i].l12;
      if (state[ch][i].l3 > l3max_t)
        l3max_t = state[ch][i].l3;

      if (state[ch][i].l12 < l12min_t)
        l12min_t = state[ch][i].l12;
      if (state[ch][i].l3 < l3min_t)
        l3min_t = state[ch][i].l3;

      if (state[ch][i].J12_2 > J12_2max_t)
        J12_2max_t = state[ch][i].J12_2;
      if (state[ch][i].j3_2 > j3_2max_t)
        j3_2max_t = state[ch][i].j3_2;

      // cout << i << "\t" << E_max_t << "\t" << E_min_t << "\t";

      if (state[ch][i].F > 0) {
        if (state[ch][i].l12 > l12max_F_t)
          l12max_F_t = state[ch][i].l12;
        if (state[ch][i].l3 > l3max_F_t)
          l3max_F_t = state[ch][i].l3;
        if (state[ch][i].l12 < l12min_F_t)
          l12min_F_t = state[ch][i].l12;
        if (state[ch][i].l3 < l3min_F_t)
          l3min_F_t = state[ch][i].l3;
        if (E_F_max_t < state[ch][i].E)
          E_F_max_t = state[ch][i].E;
        if (E_F_min_t > state[ch][i].E)
          E_F_min_t = state[ch][i].E;
        // cout << "F";
      }
      // cout << "\t";
      if (state[ch][i].I > 0) {
        if (state[ch][i].l12 > l12max_I_t)
          l12max_I_t = state[ch][i].l12;
        if (state[ch][i].l3 > l3max_I_t)
          l3max_I_t = state[ch][i].l3;
        if (state[ch][i].l12 < l12min_I_t)
          l12min_I_t = state[ch][i].l12;
        if (state[ch][i].l3 < l3min_I_t)
          l3min_I_t = state[ch][i].l3;
        if (E_I_max_t < state[ch][i].E)
          E_I_max_t = state[ch][i].E;
        if (E_I_min_t > state[ch][i].E)
          E_I_min_t = state[ch][i].E;
        // cout << "I";
      }
      // cout << endl;
    }
  }
  // cout << "*** \t " << E_F_min_t << "\t" << E_F_max_t << "\t" << E_I_min_t
  //      << "\t" << E_I_max_t << endl;
  this->l12max = l12max_t;
  this->l3max = l3max_t;
  this->l12min = l12min_t;
  this->l3min = l3min_t;
  this->lmax = std::max(l12max, l3max);

  this->l12max_F = l12max_F_t;
  this->l3max_F = l3max_F_t;
  this->l12min_F = l12min_F_t;
  this->l3min_F = l3min_F_t;

  this->l12max_I = l12max_I_t;
  this->l3max_I = l3max_I_t;
  this->l12min_I = l12min_I_t;
  this->l3min_I = l3min_I_t;

  this->n12max = n12max_t;
  this->n3max = n3max_t;
  this->nmax = std::max(n12max, n3max);
  this->n12min = n12min_t;
  this->n3min = n3min_t;
  this->nmin = std::min(n12min, n3min);
  this->J12_2max = J12_2max_t;
  this->j3_2max = j3_2max_t;
  this->all_size = all_size_t;

  this->Jaco_l123_info.clear();
  this->Jaco_l123_info.resize(2);
  // Final_state
  this->Jaco_l123_info[0].push_back(this->l12min_F);
  this->Jaco_l123_info[0].push_back(this->l12max_F);
  this->Jaco_l123_info[0].push_back(this->l3min_F);
  this->Jaco_l123_info[0].push_back(this->l3max_F);

  // Init_state
  this->Jaco_l123_info[1].push_back(this->l12min_I);
  this->Jaco_l123_info[1].push_back(this->l12max_I);
  this->Jaco_l123_info[1].push_back(this->l3min_I);
  this->Jaco_l123_info[1].push_back(this->l3max_I);

  /*cout << "++++  Jaco_l123_info information ++++" << endl;
  // cout << this->l12max_F << "\t" << this->l12max_I << endl;
  for (int i = 0; i < 2; i++) {
    for (int j = 0; j < 4; j++) {
      cout << "\t" << Jaco_l123_info[i][j];
    }
    cout << endl;
  } /**/

  this->Emin_F = E_F_min_t;
  this->Emax_F = E_F_max_t;
  this->Emin_I = E_I_min_t;
  this->Emax_I = E_I_max_t;
  this->Jaco_E_info.clear();
  this->Jaco_E_info.push_back(this->Emin_F);
  this->Jaco_E_info.push_back(this->Emax_F);
  this->Jaco_E_info.push_back(this->Emin_I);
  this->Jaco_E_info.push_back(this->Emax_I);
}

int Jacobi::FindC(int J2, int T2) {
  // Find Channel with J2 T2
  int size = channel.size();
  if (size == 0)
    return -1;
  int front = 0;
  int end = size - 1;
  int mid = size / 2 - 1;
  if (T2 == 1) {
    end = mid;
    mid = floor((end - front) / 2) + front;
  } else {
    front = mid + 1;
    mid = floor((end - front) / 2) + front;
  }
  while (front < end && channel[mid].J_2 != J2) {
    if (channel[mid].J_2 < J2)
      front = mid + 1;
    if (channel[mid].J_2 > J2)
      end = mid - 1;
    mid = (front + end) / 2;
  }
  if (channel[mid].J_2 != J2) {
    cout << "wrong at find channel in Jacobi::FindC " << endl;
    exit(0);
  } else {
    return mid;
  }
}
void Jacobi::Delete_unused_Jaco_channel(vector<int> Jaco_used_Flag) {
  int size = Jaco_used_Flag.size();
  if (size != this->state.size()) {
    cout << "Wrong happened here Jacobi::Delete_unused_Jaco_channel" << endl;
    cout << "size : " << size << "\t ch: " << this->state.size() << endl;
    exit(0);
  }
  for (int ch = 0; ch < size; ch++) {
    if (Jaco_used_Flag[ch] < 0) {
      this->state[ch].clear();
    }
  }
  this->sort();
}
//~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~end
// Jacobi~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~//

//~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~start
// Jacobi_NL~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~//
int Jacobi_NL::Tri(int &j2a, int &j2b, int &j2c) {
  // I+J>=K, I+K>=J, J+K>=I,
  // I/2+J/2+K/2 = INTEGER.
  // TRI=1, WHEN TRIADIC CONDITION IS FULFILLED, TRI=-1 OTHERWISE
  int L2 = j2a + j2b + j2c;
  if (int(L2) != (L2 / 2) * 2)
    return -1;
  L2 = L2 / 2;
  if (j2a * j2b * j2c < 0)
    return -1;
  if ((L2 - j2a) * (L2 - j2b) * (L2 - j2c) < 0)
    return -1;
  return 1;
}
void Jacobi_NL::build_NL(Jacobi &jaco, int N0_max_t, int L0_max_t) {
  N0_max = N0_max_t;
  L0_max = L0_max_t;
  J2abc_max = jaco.J2max + 2 * L0_max;

  State_Jacobi jaco_t;
  int JT_num = J2abc_max + 1;
  state_NL.resize(JT_num);
  int channel_num = -1;

  for (int T_2 = 1; T_2 <= 3; T_2 = T_2 + 2) {
    for (int J_2 = 1; J_2 <= J2abc_max; J_2 = J_2 + 2) {
      channel_num++;
      for (int i = 0; i < jaco.channel.size(); i++) {
        if (jaco.channel[i].T_2 == T_2) {
          for (int L0 = 0; L0 < L0_max; L0++) {
            int L2 = 2 * L0;
            if (Tri(L2, jaco.channel[i].J_2, J_2) == 1) {
              for (int k = 0; k < jaco.state[i].size(); k++) {
                for (int N0 = 0; N0 <= N0_max; N0++) {
                  jaco_t = jaco.state[i][k];
                  jaco_t.N_0 = N0;
                  jaco_t.L_0 = L0;
                  jaco_t.F_NL = -1;
                  jaco_t.I_NL = -1;
                  jaco_t.E_NL = jaco_t.E + 2 * N0 + L0;
                  jaco_t.Par_NL = this->phase(jaco_t.E_NL);
                  jaco_t.J2_NL = J_2;
                  jaco_t.JTNL_channel = channel_num;
                  state_NL[channel_num].push_back(jaco_t);
                }
              }
            }
          }
        }
      }
    }
  }
  channel_num++;
  this->sort();
  /*channel.resize(channel_num);
  for (int i = 0; i < channel_num; i++) {
    if (state_NL[i].size() == 0) {
      channel[i].num = 0;
      channel[i].J_2 = 0;
      channel[i].T_2 = 0;
      channel[i].E_min = 0;
      channel[i].E_max = 0;
    } else {
      channel[i].num = state_NL[i].size();
      channel[i].J_2 = state_NL[i][0].J2_NL;
      channel[i].T_2 = state_NL[i][0].T_2;
      int E_min = state_NL[i][0].E_NL;
      int E_max = state_NL[i][0].E_NL;
      for (int j = 0; j < state_NL[i].size(); j++) {
        if (E_min > state_NL[i][j].E_NL)
          E_min = state_NL[i][j].E_NL;
        if (E_max < state_NL[i][j].E_NL)
          E_max = state_NL[i][j].E_NL;
      }
      channel[i].E_min = E_min;
      channel[i].E_max = E_max;
    }
  }
  // cout<<"@@@@@@@@@@@!!!!!!!!!!!!@@@   state_NL size =
  // "<<state_NL.size()<<endl;  calulate JT2channel
  JT2channel.resize(2); // 0: T=0.5; 1:T=1.5
  for (int i = 0; i <= 1; i++) {
    int T_2 = i * 2 + 1;
    JT2channel[i].resize(J2abc_max + 1);
    for (int j = 0; j < JT2channel[i].size(); j++)
      JT2channel[i][j] = -1;
    for (int J2 = 1; J2 <= J2abc_max; J2 = J2 + 2) {
      for (int k = 0; k < channel.size(); k++) {
        if (channel[k].T_2 == T_2 && channel[k].J_2 == J2) {
          JT2channel[i][J2] = k;
          break;
        }
      }
    }
  }/**/
}

void Jacobi_NL::build_NL(Jacobi &jaco, TripleP &TP) {
  // N0_max = N0_max_t;
  // L0_max = L0_max_t;
  int J2max_TP = TP.J2_max_TP;
  // cout << "TP.J2_max_TP = " << J2max_TP << endl;
  int Ecm_Max = TP.Emax + 1;
  Emax = Ecm_Max + jaco.Emax;
  if (Emax < jaco.Emax) {
    cout << "Error in the setting of Emax in Jacobi_NL.build_NL, Emax < "
            "Jaco.Emax"
         << endl;
    exit(0);
  }
  L0_max = Ecm_Max;
  // cout << L0_max << "\t" << jaco.Emax << endl;
  // cout<<Emax<<endl;
  J2abc_max = jaco.J2max + 2 * L0_max;
  if (J2abc_max > J2max_TP)
    J2abc_max = J2max_TP;

  State_Jacobi jaco_t;
  int JT_num = J2abc_max + 1;
  state_NL.resize(JT_num);
  int channel_num = 0;
  for (int T_2 = 1; T_2 <= 3; T_2 = T_2 + 2) {
    for (int J_2 = 1; J_2 <= J2abc_max; J_2 = J_2 + 2) {
      // channel_num++;
      for (int i = 0; i < jaco.channel.size(); i++) {
        if (jaco.channel[i].T_2 == T_2) {
          for (int L0 = 0; L0 < L0_max; L0++) {
            int N0max = std::ceil((Emax - L0) / 2.0);
            int L2 = 2 * L0;
            if (Tri(L2, jaco.channel[i].J_2, J_2) == 1) {
              // cout<<"L2 = "<<L2<<"\t jaco.J_2 =
              // "<<jaco.channel[i].J_2<<"\t J_2 = "<<J_2<<endl;
              for (int k = 0; k < jaco.state[i].size(); k++) {
                for (int N0 = 0; N0 <= N0max; N0++) {
                  jaco_t = jaco.state[i][k];
                  int E = jaco_t.E + 2 * N0 + L0;
                  if (E > Emax)
                    continue;
                  //!!!! Becareful of this condition
                  // if (jaco_t.J_2 > J_2)continue;
                  jaco_t.N_0 = N0;
                  jaco_t.L_0 = L0;
                  jaco_t.F_NL = -1;
                  jaco_t.I_NL = -1;
                  jaco_t.E_NL = jaco_t.E + 2 * N0 + L0;
                  jaco_t.Par_NL = this->phase(jaco_t.E_NL);
                  jaco_t.J2_NL = J_2;
                  jaco_t.JTNL_channel = channel_num;
                  if (jaco_t.JT_channel >= jaco.channel.size()) {
                    cout << "Wrong happened @ Jacobi_NL::build_NL" << endl;
                    cout << "ch : " << channel_num << endl;
                    exit(0);
                  }
                  state_NL[channel_num].push_back(jaco_t);
                }
              }
            }
          }
        }
      }
      if (state_NL[channel_num].size() > 0)
        channel_num++;
    }
  }

  this->sort();

  /*
auto it = state_NL.begin();
while (it != state_NL.end()) {
  if ((*it).size() < 1) {
    it = state_NL.erase(it);
  } else {
    it++;
  }
}
// channel_num++;

channel.resize(state_NL.size());
for (int i = 0; i < state_NL.size(); i++) {
  if (state_NL[i].size() == 0) {
    channel[i].num = 0;
    channel[i].J_2 = 0;
    channel[i].T_2 = 0;
    channel[i].E_min = 0;
    channel[i].E_max = 0;

  } else {
    channel[i].num = state_NL[i].size();
    channel[i].J_2 = state_NL[i][0].J2_NL;
    channel[i].T_2 = state_NL[i][0].T_2;
    int E_min = state_NL[i][0].E_NL;
    int E_max = state_NL[i][0].E_NL;
    for (int j = 0; j < state_NL[i].size(); j++) {
      if (E_min > state_NL[i][j].E_NL)
        E_min = state_NL[i][j].E_NL;
      if (E_max < state_NL[i][j].E_NL)
        E_max = state_NL[i][j].E_NL;
    }
    channel[i].E_min = E_min;
    channel[i].E_max = E_max;
  }
}
// cout<<"@@@@@@@@@@@!!!!!!!!!!!!@@@   state_NL size =
// "<<state_NL.size()<<endl;  calulate JT2channel
JT2channel.resize(2); // 0: T=0.5; 1:T=1.5
for (int i = 0; i <= 1; i++) {
  int T_2 = i * 2 + 1;
  JT2channel[i].resize(J2abc_max + 1);
  for (int j = 0; j < JT2channel[i].size(); j++)
    JT2channel[i][j] = -1;
  for (int J2 = 1; J2 <= J2abc_max; J2 = J2 + 2) {
    for (int k = 0; k < channel.size(); k++) {
      if (channel[k].T_2 == T_2 && channel[k].J_2 == J2) {
        JT2channel[i][J2] = k;
        break;
      }
    }
  }
}/**/
}

void Jacobi_NL::sort() {
  auto it = state_NL.begin();
  while (it != state_NL.end()) {
    if ((*it).size() < 1) {
      it = state_NL.erase(it);
    } else {
      it++;
    }
  }
  int ch_size = state_NL.size();
  // sort by E
  for (int i = 0; i < ch_size; i++)
    stable_sort(this->state_NL[i].begin(), this->state_NL[i].end(),
                [](State_Jacobi a, State_Jacobi b) {
                  return a.E_NL * 100000 + a.N_0 * 100 + a.L_0 +
                             0.01 * a.JT_channel + 0.0001 * a.index <=
                         b.E_NL * 100000 + b.N_0 * 100 + b.L_0 +
                             0.01 * b.JT_channel + 0.0001 * b.index;
                });
  for (int ch = 0; ch < ch_size; ch++) {
    int index_NL_t = 0;
    for (int i = 0; i < state_NL[ch].size(); i++) {
      state_NL[ch][i].index_NL = index_NL_t;
      state_NL[ch][i].JTNL_channel = ch;
      index_NL_t++;
    }
  }
  channel.clear();
  channel.resize(state_NL.size());
  for (int i = 0; i < state_NL.size(); i++) {
    if (state_NL[i].size() == 0) {
      channel[i].num = 0;
      channel[i].J_2 = 0;
      channel[i].T_2 = 0;
      channel[i].E_min = 0;
      channel[i].E_max = 0;

    } else {
      channel[i].num = state_NL[i].size();
      channel[i].J_2 = state_NL[i][0].J2_NL;
      channel[i].T_2 = state_NL[i][0].T_2;
      int E_min = state_NL[i][0].E_NL;
      int E_max = state_NL[i][0].E_NL;
      for (int j = 0; j < state_NL[i].size(); j++) {
        if (E_min > state_NL[i][j].E_NL)
          E_min = state_NL[i][j].E_NL;
        if (E_max < state_NL[i][j].E_NL)
          E_max = state_NL[i][j].E_NL;
      }
      channel[i].E_min = E_min;
      channel[i].E_max = E_max;
    }
  }

  int all_size_t = 0;
  for (int i = 0; i < state_NL.size(); i++) {
    all_size_t += state_NL[i].size();
    int index_NL = 0;
    int E_min = channel[i].E_min;
    int E_max = channel[i].E_max;
    channel[i].indexE.resize(E_max + 1);
    channel[i].indexE_beg.resize(E_max + 1);
    for (int k = 0; k < state_NL[i].size(); k++) {
      state_NL[i][k].index_NL = index_NL;
      channel[i].indexE[state_NL[i][k].E_NL] = index_NL;
      index_NL++;
    }
    for (int E = E_min; E <= E_max; E++) {
      int indexE_beg = 0;
      for (int k = 0; k < state_NL[i].size(); k++) {
        if (state_NL[i][k].E_NL == E) {
          indexE_beg = k;
          break;
        }
      }
      channel[i].indexE_beg[E] = indexE_beg;
    }
  }
  this->all_size = all_size_t;
  // find channel_size_T1
  int size_T1 = 0;
  for (int i = 0; i < channel.size(); i++) {
    if (channel[i].T_2 == 1) {
      size_T1++;
    }
  }
  this->channel_size_T1 = size_T1;
  JT2channel.resize(2); // 0: T=0.5; 1:T=1.5
  for (int i = 0; i <= 1; i++) {
    int T_2 = i * 2 + 1;
    JT2channel[i].resize(J2abc_max + 1);
    for (int j = 0; j < JT2channel[i].size(); j++)
      JT2channel[i][j] = -1;
    for (int J2 = 1; J2 <= J2abc_max; J2 = J2 + 2) {
      for (int k = 0; k < channel.size(); k++) {
        if (channel[k].T_2 == T_2 && channel[k].J_2 == J2) {
          JT2channel[i][J2] = k;
          break;
        }
      }
    }
  }
}
int Jacobi_NL::FindC(int J2, int T2) {
  // Find Channel with J2 T2
  int size = channel.size();
  if (size == 0)
    return -1;
  int front = 0;
  int end = size - 1;
  int mid = this->channel_size_T1 - 1;
  if (T2 == 1) {
    end = mid;
    mid = floor((end - front) / 2) + front;
  } else {
    front = mid + 1;
    mid = floor((end - front) / 2) + front;
  }
  while (front < end && channel[mid].J_2 != J2) {
    if (channel[mid].J_2 < J2)
      front = mid + 1;
    if (channel[mid].J_2 > J2)
      end = mid - 1;
    mid = (front + end) / 2;
  }
  if (channel[mid].J_2 != J2) {
    cout << "wrong at find channel in Jacobi_NL::FindC " << endl;
    cout << " J2 = " << J2 << " T2 = " << T2 << endl;
    exit(0);
  } else {
    return mid;
  }
}
void Jacobi_NL::select_FI_T(TripleP &TP) {
  for (int ch = 0; ch < state_NL.size(); ch++) {
    int J2 = this->channel[ch].J_2;
    int T2 = this->channel[ch].T_2;
    // cout << "J2 = " << J2 << "\t T2 = ";
    int ch_TP = TP.FindC(J2, T2);
    // cout << "ch_TP = " << ch_TP << endl;
    if (ch_TP < 0)
      continue;
    vector<int> E_F_vec, E_I_vec;
    int E_f_t = -1;
    int E_i_t = -1;
    for (int k = 0; k < TP.TPconf_E.size(); k++) {
      for (int j = 0; j < TP.TPconf_E[k][ch_TP].size(); j++) {
        if (E_f_t < TP.TPconf_E[k][ch_TP][j].E_F) {
          E_f_t = TP.TPconf_E[k][ch_TP][j].E_F;
          E_F_vec.push_back(E_f_t);
        }
        if (E_i_t < TP.TPconf_E[k][ch_TP][j].E_I) {
          E_i_t = TP.TPconf_E[k][ch_TP][j].E_I;
          E_I_vec.push_back(E_i_t);
        }
      }
    }
    for (int i = 0; i < state_NL[ch].size(); i++) {
      int E = state_NL[ch][i].E_NL;
      auto it_f = find(E_F_vec.begin(), E_F_vec.end(), E);
      auto it_i = find(E_I_vec.begin(), E_I_vec.end(), E);
      if (it_f != E_F_vec.end()) {
        state_NL[ch][i].F_NL = 1;
      }
      if (it_i != E_I_vec.end()) {
        state_NL[ch][i].I_NL = 1;
      }
    }
  }
  // cout << "Begin to delete Jacobi_NL" << endl;
  // !!! Delete the state
  /*
   for (int ch = 0; ch < state_NL.size(); ch++) {
     cout << "ch = " << ch << "\t size : " << state_NL[ch].size();
     auto it = state_NL[ch].begin();
     while (it != state_NL[ch].end()) {
       if (it->F_NL < 0 && it->I_NL < 0) {
         it = state_NL[ch].erase(it);
       } else {
         it++;
       }
     }
     cout << "\t delete size : " << state_NL[ch].size() << endl;
   } /**/
  // cout << "End to delete Jacobi_NL" << endl;
}
void Jacobi_NL::Delet_Sort_T() {
  vector<vector<State_Jacobi>> state_NL_T; // state[channel][i]
  state_NL_T = state_NL;
  state_NL.clear();
  state_NL.resize(state_NL_T.size());
  for (int ch = 0; ch < state_NL_T.size(); ch++) {
    for (int i = 0; i < state_NL_T[ch].size(); i++) {
      if (state_NL_T[ch][i].F_NL < 0 && state_NL_T[ch][i].I_NL < 0)
        continue;
      state_NL[ch].push_back(state_NL_T[ch][i]);
    }
  }
  state_NL_T.clear();
}
void Jacobi_NL::Delete_unused_Jaco_channel(Jacobi &jaco,
                                           vector<int> &Jaco_used_Flag) {
  // vector<int> Jaco_used_Flag;
  int jaco_ch_size = jaco.state.size();
  // cout << "\t !! ## !! jaco.state.size() : " << jaco.state.size() << endl;
  int jacoNL_ch_size = this->state_NL.size();
  Jaco_used_Flag.resize(jaco_ch_size);
  for (int i = 0; i < jaco_ch_size; i++) {
    Jaco_used_Flag[i] = -1;
  }
  for (int ch = 0; ch < jacoNL_ch_size; ch++) {
    for (int i = 0; i < state_NL[ch].size(); i++) {
      int ch_jaco = state_NL[ch][i].JT_channel;
      // cout<<"ch_jaco = "<<ch_jaco<<endl;
      if (ch_jaco >= jaco_ch_size) {
        cout << "Wrong happened Jacobi_NL::Delete_unused_Jaco_channel" << endl;
        cout << "ch_NL : " << ch << "\t i : " << i
             << "\t index : " << state_NL[ch][i].index << endl;
        cout << "JT_channel : " << ch_jaco << "\t ch_jaco : " << jaco_ch_size
             << endl;
        exit(0);
      }
      Jaco_used_Flag[ch_jaco] = 1;
    }
  }
  // cout << "\t Jaco_used_Flag.size() = " << Jaco_used_Flag.size() << endl;
}
void Jacobi_NL::print() {
  // cout<<"total number of basis_JT is : "<<endl;
  // cout<<num<<endl;
  for (int i = 0; i < channel.size(); i++) {
    cout << "===JT channel : " << i << " ====" << endl;
    cout << "num: \t"
         << "2*T \t"
         << "2*T12 \t"
         << "2*J \t"
         << "2*J12 \t"
         << "2*j3 \t"
         << "2*S12 \t"
         << "l12 \t"
         << "l3 \t"
         << "n12 \t"
         << "n3 \t"
         << "E \t"
         << "index \t"
         << "JT_channel \t"
         << "flag \t"
         << "J2_NL \t"
         << "L0 \t"
         << "E_NL \t"
         << "JTNL_channel \t"
         << "index_NL" << endl;
    for (int j = 0; j < state_NL[i].size(); j++) {
      cout << j + 1 << "\t";
      cout << state_NL[i][j].T_2 << "\t" << state_NL[i][j].T12_2 << "\t"
           << state_NL[i][j].J_2 << "\t" << state_NL[i][j].J12_2 << "\t";
      cout << state_NL[i][j].j3_2 << "\t";
      cout << state_NL[i][j].S12_2 << "\t" << state_NL[i][j].l12 << "\t"
           << state_NL[i][j].l3 << "\t" << state_NL[i][j].n12 << "\t";
      cout << state_NL[i][j].n3 << "\t" << state_NL[i][j].E << "\t"
           << state_NL[i][j].index << "\t";
      cout << state_NL[i][j].JT_channel << "\t\t" << state_NL[i][j].flag << "\t"
           << state_NL[i][j].J2_NL << "\t";
      cout << state_NL[i][j].L_0 << "\t" << state_NL[i][j].E_NL << " \t"
           << state_NL[i][j].JTNL_channel << " \t" << state_NL[i][j].index_NL
           << endl;
    }
  }
}

void Jacobi_NL::print(int channel_num) {
  // cout<<"total number of basis_JT is : "<<endl;
  // cout<<num<<endl;
  int i = channel_num;
  // for(int i=0;i<channel.size();i++){
  cout << "===JT channel : " << i << " ====" << endl;
  cout << "num: \t"
       << "2*T \t"
       << "2*T12 \t"
       << "2*J \t"
       << "2*J12 \t"
       << "2*j3 \t"
       << "2*S12 \t"
       << "l12 \t"
       << "l3 \t"
       << "n12 \t"
       << "n3 \t"
       << "E \t"
       << "index \t"
       << "JT_channel \t"
       << "flag \t"
       << "J2_NL \t"
       << "L0 \t"
       << "E_NL \t"
       << "JTNL_channel \t"
       << "index_NL" << endl;
  for (int j = 0; j < state_NL[i].size(); j++) {
    cout << j + 1 << "\t";
    cout << state_NL[i][j].T_2 << "\t" << state_NL[i][j].T12_2 << "\t"
         << state_NL[i][j].J_2 << "\t" << state_NL[i][j].J12_2 << "\t";
    cout << state_NL[i][j].j3_2 << "\t";
    cout << state_NL[i][j].S12_2 << "\t" << state_NL[i][j].l12 << "\t"
         << state_NL[i][j].l3 << "\t" << state_NL[i][j].n12 << "\t";
    cout << state_NL[i][j].n3 << "\t" << state_NL[i][j].E << "\t"
         << state_NL[i][j].index << "\t";
    cout << state_NL[i][j].JT_channel << "\t\t" << state_NL[i][j].flag << "\t"
         << state_NL[i][j].J2_NL << "\t";
    cout << state_NL[i][j].L_0 << "\t" << state_NL[i][j].E_NL << " \t"
         << state_NL[i][j].JTNL_channel << " \t" << state_NL[i][j].index_NL
         << endl;
  }
  //}
}

void Jacobi_NL::print(std::string filename) {
  std::ofstream file;
  file.open(filename);
  // file<<"Jacobi_NL basis with L0 max = "<<L0_max<<", \t N0 max =
  // "<<N0_max<<endl;
  file << "Energy cut off is : " << Emax << " ( hbar-omega )" << endl;
  file << "total number of basis is : " << endl;
  int total_num = 0;
  for (int i = 0; i < state_NL.size(); i++) {
    total_num += state_NL[i].size();
  }
  file << total_num << endl;
  // file<<"total number of basis is : "<<endl;
  // file<<this->num<<endl;
  file << state_NL.size() << "\t JT channel_NL num" << endl;
  for (int i = 0; i < state_NL.size(); i++) {
    file << state_NL[i].size() << "  ";
  }
  file << endl;
  file << "num: \t"
       << "2*T \t"
       << "2*T12 \t"
       << "2*J123 \t"
       << "2*J12 \t"
       << "2*j3 \t"
       << "2*S12 \t"
       << "l12 \t"
       << "l3 \t"
       << "n12 \t"
       << "n3 \t"
       << "E \t"
       << "index \t"
       << "Jaco_channel \t"
       << "flag \t"
       << "2*J_NL \t"
       << "N0 \t"
       << "L0 \t"
       << "E_NL \t"
       << "index_NL\t F\t I" << endl;
  for (int i = 0; i < state_NL.size(); i++) {
    file << "channel_NL No. : " << i << endl;
    for (int j = 0; j < state_NL[i].size(); j++) {
      file << j + 1 << "\t";
      file << state_NL[i][j].T_2 << "\t" << state_NL[i][j].T12_2 << "\t"
           << state_NL[i][j].J_2 << "\t" << state_NL[i][j].J12_2 << "\t";
      file << state_NL[i][j].j3_2 << "\t";
      file << state_NL[i][j].S12_2 << "\t" << state_NL[i][j].l12 << "\t"
           << state_NL[i][j].l3 << "\t" << state_NL[i][j].n12 << "\t";
      file << state_NL[i][j].n3 << "\t" << state_NL[i][j].E << "\t"
           << state_NL[i][j].index << "\t";
      file << state_NL[i][j].JT_channel << "\t\t" << state_NL[i][j].flag << "\t"
           << state_NL[i][j].J2_NL << "\t";
      file << state_NL[i][j].N_0 << "\t" << state_NL[i][j].L_0 << "\t"
           << state_NL[i][j].E_NL << "\t"
           << state_NL[i][j].index_NL; // << "\t ";
      /*if (state_NL[i][j].F_NL > 0)
        file << "\t *";
      file << "\t ";
      if (state_NL[i][j].I_NL > 0)
        file << "*";/**/
      file << endl;
    }
  }
  file.close();
}

void Jacobi_NL::print_short(std::string filename) {
  std::ofstream file;
  file.open(filename);
  // file<<"Jacobi_NL basis with L0 max = "<<L0_max<<", \t N0 max =
  // "<<N0_max<<endl;
  file << "Energy cut off is : " << Emax << " ( hbar-omega )" << endl;
  file << "total number of basis is : " << endl;
  int total_num = 0;
  for (int i = 0; i < state_NL.size(); i++) {
    total_num += state_NL[i].size();
  }
  file << total_num << endl;
  // file<<"total number of basis is : "<<endl;
  // file<<this->num<<endl;
  file << state_NL.size() << "\t JT channel_NL num" << endl;
  for (int i = 0; i < state_NL.size(); i++) {
    file << state_NL[i].size() << "  ";
  }
  file << endl;
  file << "num: \t"
       << "2*T \t"
       << "2*T12 \t"
       << "2*J123 \t"
       << "2*J12 \t"
       << "2*j3 \t"
       << "2*S12 \t"
       << "l12 \t"
       << "l3 \t"
       << "n12 \t"
       << "n3 \t"
       << "2*J_NL \t"
       << "N0 \t"
       << "L0 \t"
       << "E_NL \t"
       << "index_NL" << endl;
  for (int i = 0; i < state_NL.size(); i++) {
    file << "channel_NL No. : " << i << endl;
    for (int j = 0; j < state_NL[i].size(); j++) {
      file << j + 1 << "\t";
      file << state_NL[i][j].T_2 << "\t" << state_NL[i][j].T12_2 << "\t"
           << state_NL[i][j].J_2 << "\t" << state_NL[i][j].J12_2 << "\t";
      file << state_NL[i][j].j3_2 << "\t";
      file << state_NL[i][j].S12_2 << "\t" << state_NL[i][j].l12 << "\t"
           << state_NL[i][j].l3 << "\t" << state_NL[i][j].n12 << "\t";
      file << state_NL[i][j].n3 << "\t";
      file << state_NL[i][j].J2_NL << "\t";
      file << state_NL[i][j].N_0 << "\t" << state_NL[i][j].L_0 << "\t"
           << state_NL[i][j].E_NL << "\t" << state_NL[i][j].index_NL << endl;
    }
  }
  file.close();
}

void Jacobi_NL::print_channel() {
  int size = this->channel.size();
  cout << " Jacobi_NL channel Total Num: " << size << endl;
  cout << "index \t"
       << "num \t"
       << "J2 \t"
       << "T_2 \t"
       << "E_min \t"
       << "E_max \t" << endl;
  for (int i = 0; i < size; i++) {
    cout << i << "\t" << channel[i].num << "\t" << channel[i].J_2 << "\t"
         << channel[i].T_2 << "\t" << channel[i].E_min << "\t"
         << channel[i].E_max << endl;
  }
}
void Jacobi_NL::print_channel(std::string filename) {
  std::ofstream file;
  file.open(filename);
  int size = this->channel.size();
  file << " Jacobi_NL channel Total Num: " << size << endl;
  file << "index \t"
       << "num \t"
       << "J2 \t"
       << "T_2 \t"
       << "E_min \t"
       << "E_max \t" << endl;
  for (int i = 0; i < size; i++) {
    file << i << "\t" << channel[i].num << "\t" << channel[i].J_2 << "\t"
         << channel[i].T_2 << "\t" << channel[i].E_min << "\t"
         << channel[i].E_max << endl;
  }
  file << "E_min and E_max index:" << endl;
  file << "\t \t   "
       << "E : 0\t1\t2 \t3 \t4 \t5 ....." << endl;
  for (int i = 0; i < size; i++) {
    file << "channel : " << i << endl;
    file << "\t E_min index : ";
    for (int j = 0; j < channel[i].indexE_beg.size(); j++) {
      file << channel[i].indexE_beg[j] << "\t";
    }
    file << endl;
    file << "\t E_max index : ";
    for (int j = 0; j < channel[i].indexE.size(); j++) {
      file << channel[i].indexE[j] << "\t";
    }
    file << endl;
  }
  file.close();
}
//~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~end
// Jacobi_NL~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~//
//~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~start Single
// Particle~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~//
void SingleP::build_T(std::string filename) {
  std::ifstream file;
  file.open(filename.data());
  if (!file) {
    std::cout << "wrong at SP file path" << std::endl;
    exit(0);
  }

  // state.clear();

  char u1[256];   // unuse
  std::string u2; // unuse
  double u3;
  int n_max_t, l_max_t;
  int num;

  // cout << filename << endl;
  file.getline(u1, 256, '\n');
  file >> u2 >> u2 >> n_max_t >> u2 >> u2 >> l_max_t;
  file.getline(u1, 256, '\n');
  file >> num;
  file.getline(u1, 256, '\n');
  file.getline(u1, 256, '\n');

  // n_max = n_max_t;
  // l_max = l_max_t;
  State_SP state_t;
  // int index=0;
  // cout << "n_max_t = " << n_max << endl;
  // cout << "l_max_t = " << l_max << endl;
  // cout << "Read SP Num :" << num << endl;
  for (int i = 0; i < num; i++) {
    file >> state_t.tz_2 >> state_t.n >> state_t.l >> state_t.j_2 >>
        state_t.E >> state_t.index;
    // state_t.index=index;
    file.getline(u1, 256, '\n');
    // index++;
    this->state.push_back(state_t);
  }
  file.close();
  n_max_t = -1;
  l_max_t = -1;
  for (int i = 0; i < num; i++) {
    if (state[i].n > n_max_t)
      n_max_t = state[i].n;
    if (state[i].l > l_max_t)
      l_max_t = state[i].l;
  }
  n_max = n_max_t;
  l_max = l_max_t;
}

void SingleP::build_T(vector<vector<int>> sp_vec) {
  cout<<"   Start build_T(vector<vector<int>> sp_vec) "<<endl;
  // state.clear();

  char u1[256];   // unuse
  std::string u2; // unuse
  double u3;
  int n_max_t, l_max_t;
  int num = sp_vec.size();

  // cout << filename << endl;
  // file.getline(u1, 256, '\n');
  // file >> u2 >> u2 >> n_max_t >> u2 >> u2 >> l_max_t;
  // file.getline(u1, 256, '\n');
  // file >> num;
  // file.getline(u1, 256, '\n');
  // file.getline(u1, 256, '\n');

  // n_max = n_max_t;
  // l_max = l_max_t;
  State_SP state_t;
  // int index=0;
  // cout << "n_max_t = " << n_max << endl;
  // cout << "l_max_t = " << l_max << endl;
  // cout << "Read SP Num :" << num << endl;
  for (int i = 0; i < num; i++) {

    state_t.tz_2 = sp_vec[i][0];
    state_t.n = sp_vec[i][1];
    state_t.l = sp_vec[i][2];
    state_t.j_2 = sp_vec[i][3];
    state_t.E = sp_vec[i][4];
    state_t.index = sp_vec[i][5];
    // index++;
    this->state.push_back(state_t);
  }
  // file.close();
  n_max_t = -1;
  l_max_t = -1;
  for (int i = 0; i < num; i++) {
    if (state[i].n > n_max_t)
      n_max_t = state[i].n;
    if (state[i].l > l_max_t)
      l_max_t = state[i].l;
  }
  this->n_max = n_max_t;
  this->l_max = l_max_t;
}

void SingleP::build_T_E(std::string filename) {
  int Emax = -1;
  std::ifstream file;
  file.open(filename.data());
  if (!file) {
    std::cout << "wrong at SP file path E" << std::endl;
    exit(0);
  }
  char u1[256];   // unuse
  std::string u2; // unuse
  file >> Emax;
  file.getline(u1, 256, '\n');
  file.getline(u1, 256, '\n');

  int N_exch = 0;
  vector<int> index_A, index_B;
  if (!file.eof()) {
    file >> N_exch;
    cout << "\t  -- N_exch = " << N_exch << endl;
    file.getline(u1, 256, '\n');
    index_A.resize(N_exch);
    index_B.resize(N_exch);
    for (int i = 0; i < N_exch; i++) {
      int a, b;
      file >> a >> b;
      index_A[i] = a;
      index_B[i] = b;
      file.getline(u1, 256, '\n');
      if (file.eof())
        break;
    }
  }
  file.close();
  cout << " SingleP::build_T_E E = " << Emax << endl;
  // n_max = n_max_t;
  // l_max = l_max_t;
  int n_max_t = -1;
  int l_max_t = -1;

  State_SP state_t;
  // int index=0;
  for (int E = 0; E <= Emax; E++) {
    for (int n = 0; n <= (Emax / 2 + 2); n++) {
      for (int l = 0; l <= Emax; l++) {
        if ((2 * n + l) != E)
          continue;
        int j2_min = std::abs(2 * l - 1);
        int j2_max = std::abs(2 * l + 1);
        for (int j2 = j2_min; j2 <= j2_max; j2 = j2 + 2) {

          state_t.n = n;
          state_t.l = l;
          state_t.j_2 = j2;
          state_t.tz_2 = 0;
          state_t.E = 2 * n + l;
          // state_t.index=index;
          // index++;
          this->state.push_back(state_t);
          if (n > n_max_t)
            n_max_t = n;
          if (l > l_max_t)
            l_max_t = l;
        }
      }
    }
  }
  this->n_max = n_max_t;
  this->l_max = l_max_t;
  stable_sort(state.begin(), state.end(), [](State_SP a, State_SP b) {
    return (10000 * a.E + 100 * a.l - a.j_2 + 0.1 * a.tz_2) <=
           (10000 * b.E + 100 * b.l - b.j_2 + 0.1 * b.tz_2);
  });
  int index = 0;
  for (int i = 0; i < state.size(); i++) {
    state[i].index = index;
    index++;
  }
  if (N_exch > 0) {
    for (int i = 0; i < N_exch; i++) {
      int A = index_A[i];
      int B = index_B[i];
      state[A].index = B;
      state[B].index = A;
      cout << "\t --- A:" << A << "\t B:" << B << endl;
    }
    stable_sort(state.begin(), state.end(),
                [](State_SP a, State_SP b) { return (a.index) <= (b.index); });
  }
}

void SingleP::build_T_E(int Emax) {
  cout << " SingleP::build_T_E E = " << Emax << endl;
  // n_max = n_max_t;
  // l_max = l_max_t;
  int n_max_t = -1;
  int l_max_t = -1;

  State_SP state_t;
  // int index=0;
  for (int E = 0; E <= Emax; E++) {
    for (int n = 0; n <= (Emax / 2 + 2); n++) {
      for (int l = 0; l <= Emax; l++) {
        if ((2 * n + l) != E)
          continue;
        int j2_min = std::abs(2 * l - 1);
        int j2_max = std::abs(2 * l + 1);
        for (int j2 = j2_min; j2 <= j2_max; j2 = j2 + 2) {

          state_t.n = n;
          state_t.l = l;
          state_t.j_2 = j2;
          state_t.tz_2 = 0;
          state_t.E = 2 * n + l;
          // state_t.index=index;
          // index++;
          this->state.push_back(state_t);
          if (n > n_max_t)
            n_max_t = n;
          if (l > l_max_t)
            l_max_t = l;
        }
      }
    }
  }
  this->n_max = n_max_t;
  this->l_max = l_max_t;
  stable_sort(state.begin(), state.end(), [](State_SP a, State_SP b) {
    return (10000 * a.E + 100 * a.l - a.j_2 + 0.1 * a.tz_2) <=
           (10000 * b.E + 100 * b.l - b.j_2 + 0.1 * b.tz_2);
  });
  int index = 0;
  for (int i = 0; i < state.size(); i++) {
    state[i].index = index;
    index++;
  }
}

void SingleP::build_Tz(std::string filename) {
  std::ifstream file;
  file.open(filename.data());
  if (!file) {
    std::cout << "wrong at SP file path" << std::endl;
    exit(0);
  }

  // state.clear();

  char u1[256];   // unuse
  std::string u2; // unuse
  double u3;
  int n_max_t, l_max_t;
  int num;

  cout << filename << endl;
  file.getline(u1, 256, '\n');
  file >> u2 >> u2 >> n_max_t >> u2 >> u2 >> l_max_t;
  file.getline(u1, 256, '\n');
  file >> num;
  file.getline(u1, 256, '\n');
  file.getline(u1, 256, '\n');

  n_max = n_max_t;
  l_max = l_max_t;
  State_SP state_t;
  // int index=0;
  // cout << "n_max_t = " << n_max << endl;
  // cout << "l_max_t = " << l_max << endl;
  // cout << "Read SP Num :" << num << endl;
  this->state_Tz.clear();
  for (int i = 0; i < num; i++) {
    file >> state_t.tz_2 >> state_t.n >> state_t.l >> state_t.j_2 >>
        state_t.E >> state_t.index;
    // state_t.index=index;
    /*cout << state_t.tz_2 << "\t" << state_t.n << "\t" << state_t.l << "\t"
         << state_t.j_2 << "\t" << state_t.E << "\t" << state_t.index <<
       endl;/**/
    file.getline(u1, 256, '\n');
    // index++;
    this->state_Tz.push_back(state_t);
  }
  file.close();
}
void SingleP::build_Tz(int n_max_t, int l_max_t) {
  n_max = n_max_t;
  l_max = l_max_t;
  State_SP state_t;
  // int index=0;
  for (int tz = -1; tz <= 1; tz = tz + 2) {
    for (int n = 0; n <= n_max; n++) {
      for (int l = 0; l <= l_max; l++) {
        int j2_min = std::abs(2 * l - 1);
        int j2_max = std::abs(2 * l + 1);
        for (int j2 = j2_min; j2 <= j2_max; j2 = j2 + 2) {
          state_t.n = n;
          state_t.l = l;
          state_t.j_2 = j2;
          state_t.tz_2 = tz;
          state_t.E = 2 * n + l;
          // state_t.index=index;
          // index++;
          this->state_Tz.push_back(state_t);
        }
      }
    }
  }
  stable_sort(state_Tz.begin(), state_Tz.end(), [](State_SP a, State_SP b) {
    return (10000 * a.E + 100 * a.l - a.j_2 + 0.1 * a.tz_2) <=
           (10000 * b.E + 100 * b.l - b.j_2 + 0.1 * b.tz_2);
  });
  int index = 0;
  for (int i = 0; i < state_Tz.size(); i++) {
    state_Tz[i].index = index;
    index++;
  }
}
void SingleP::build_Tz(int n_max_t, int l_max_t, int l_bad, int j2_bad) {
  n_max = n_max_t;
  l_max = l_max_t;
  State_SP state_t;
  // int index=0;
  for (int tz = -1; tz <= 1; tz = tz + 2) {
    for (int n = 0; n <= n_max; n++) {
      for (int l = 0; l <= l_max; l++) {
        int j2_min = std::abs(2 * l - 1);
        int j2_max = std::abs(2 * l + 1);
        for (int j2 = j2_min; j2 <= j2_max; j2 = j2 + 2) {
          if (l == l_bad && j2 == j2_bad)
            continue;
          state_t.n = n;
          state_t.l = l;
          state_t.j_2 = j2;
          state_t.tz_2 = tz;
          state_t.E = 2 * n + l;
          // state_t.index=index;
          // index++;
          this->state_Tz.push_back(state_t);
        }
      }
    }
  }
  stable_sort(state_Tz.begin(), state_Tz.end(), [](State_SP a, State_SP b) {
    return (10000 * a.E + 100 * a.l - a.j_2 + 0.1 * a.tz_2) <=
           (10000 * b.E + 100 * b.l - b.j_2 + 0.1 * b.tz_2);
  });
  int index = 0;
  for (int i = 0; i < state_Tz.size(); i++) {
    state_Tz[i].index = index;
    index++;
  }
}
void SingleP::build_Tz() {
  // cout << "Becareful !!! here use SingleP::build build state to build "
  //        "state_Tz "
  //     << endl;
  if (state.size() < 1) {
    cout << "build_T first !" << endl;
    exit(0);
  }
  this->state_Tz.clear();
  for (int i = 0; i < state.size(); i++) {
    State_SP state_t;
    state_t.n = state[i].n;
    state_t.l = state[i].l;
    state_t.j_2 = state[i].j_2;
    state_t.E = state[i].E;
    state_t.tz_2 = -1;
    this->state_Tz.push_back(state_t);
    state_t.tz_2 = 1;
    this->state_Tz.push_back(state_t);
  }
  int index = 0;
  for (int i = 0; i < state_Tz.size(); i++) {
    state_Tz[i].index = index;
    index++;
  }
}

void SingleP::build_TzM(int n_max_t, int l_max_t, int l_bad, int j2_bad) {
  n_max = n_max_t;
  l_max = l_max_t;
  State_SPM state_t;
  // int index=0;
  for (int tz = -1; tz <= 1; tz = tz + 2) {
    for (int n = 0; n <= n_max; n++) {
      for (int l = 0; l <= l_max; l++) {
        int j2_min = std::abs(2 * l - 1);
        int j2_max = std::abs(2 * l + 1);
        for (int j2 = j2_min; j2 <= j2_max; j2 = j2 + 2) {
          if (l == l_bad && j2 == j2_bad)
            continue;
          for (int m_2 = -1 * j2; m_2 <= j2; m_2 = m_2 + 2) {
            state_t.n = n;
            state_t.l = l;
            state_t.j_2 = j2;
            state_t.m_2 = m_2;
            state_t.tz_2 = tz;
            state_t.E = 2 * n + l;
            // state_t.index=index;
            // index++;
            this->state_TzM.push_back(state_t);
          }
        }
      }
    }
  }
  stable_sort(state_TzM.begin(), state_TzM.end(), [](State_SPM a, State_SPM b) {
    return (100 * a.E + 10 * a.l - a.j_2 + 1000 * a.tz_2 + 0.01 * a.m_2) <=
           (100 * b.E + 10 * b.l - b.j_2 + 1000 * b.tz_2 + 0.01 * b.m_2);
  });
  int index = 0;
  for (int i = 0; i < state_TzM.size(); i++) {
    state_TzM[i].index = index;
    index++;
  }
  int size = state_TzM.size();
  index_M2J.resize(size);
  for (int i = 0; i < size; i++) {
    int n = state_TzM[i].n;
    int l = state_TzM[i].l;
    int j_2 = state_TzM[i].j_2;
    int m_2 = state_TzM[i].m_2;
    int tz_2 = state_TzM[i].tz_2;
    int E = state_TzM[i].E;
    int index = -1;
    for (int j = 0; j < state_Tz.size(); j++) {
      if (state_Tz[j].n != n || state_Tz[j].l != l)
        continue;
      if (state_Tz[j].j_2 != j_2)
        continue;
      if (state_Tz[j].tz_2 != tz_2 || state_Tz[j].E != E)
        continue;
      index = j;
      break;
    }
    index_M2J[i] = index;
    // cout<<i<<"\t"<<index<<endl;
  }
}

void SingleP::build_T(int n_max_t, int l_max_t, int l_bad, int j2_bad) {
  n_max = n_max_t;
  l_max = l_max_t;
  State_SP state_t;
  // int index=0;
  for (int tz = 0; tz <= 0; tz = tz + 2) {
    for (int n = 0; n <= n_max; n++) {
      for (int l = 0; l <= l_max; l++) {
        int j2_min = std::abs(2 * l - 1);
        int j2_max = std::abs(2 * l + 1);
        for (int j2 = j2_min; j2 <= j2_max; j2 = j2 + 2) {
          if (l == l_bad && j2 == j2_bad)
            continue;
          state_t.n = n;
          state_t.l = l;
          state_t.j_2 = j2;
          state_t.tz_2 = tz;
          state_t.E = 2 * n + l;
          // state_t.index=index;
          // index++;
          this->state.push_back(state_t);
        }
      }
    }
  }
  stable_sort(state.begin(), state.end(), [](State_SP a, State_SP b) {
    return (10000 * a.E + 100 * a.l - a.j_2 + 0.1 * a.tz_2) <=
           (10000 * b.E + 100 * b.l - b.j_2 + 0.1 * b.tz_2);
  });
  int index = 0;
  for (int i = 0; i < state.size(); i++) {
    state[i].index = index;
    index++;
  }
}

void SingleP::print() {
  cout << "single particle state" << endl;
  cout << "n_max :" << n_max << "\t l_max" << l_max << endl;
  cout << "tz \t"
       << "n \t"
       << "l \t"
       << "j2 \t"
       << "E \t"
       << "index" << endl;
  for (int i = 0; i < state.size(); i++) {
    cout << state[i].tz_2 << "\t";
    cout << state[i].n << "\t";
    cout << state[i].l << "\t";
    cout << state[i].j_2 << "\t";
    cout << state[i].E << "\t";
    cout << state[i].index << "\t";
    cout << endl;
  }
}
void SingleP::print(int a_min, int a_max, int b_min, int b_max, int c_min,
                    int c_max) {
  cout << "single particle state" << endl;
  cout << "n_max :" << n_max << "\t l_max" << l_max << endl;
  cout << "tz \t"
       << "n \t"
       << "l \t"
       << "j2 \t"
       << "E \t"
       << "index \t"
       << "a \t"
       << "b \t"
       << "c \t" << endl;
  for (int i = 0; i < state.size(); i++) {
    cout << state[i].tz_2 << "\t";
    cout << state[i].n << "\t";
    cout << state[i].l << "\t";
    cout << state[i].j_2 << "\t";
    cout << state[i].E << "\t";
    cout << state[i].index << "\t";
    if (i >= a_min && i <= a_max) {
      cout << "*\t";
    } else {
      cout << " \t";
    }
    if (i >= b_min && i <= b_max) {
      cout << "*\t";
    } else {
      cout << " \t";
    }
    if (i >= c_min && i <= c_max) {
      cout << "*\t";
    } else {
      cout << " \t";
    }
    cout << endl;
  }
}
void SingleP::print(vector<vector<int>> TP_allow) {
  cout << "single particle state" << endl;
  cout << "n_max :" << n_max << "\t l_max" << l_max << endl;
  cout << "tz \t"
       << "n \t"
       << "l \t"
       << "j2 \t"
       << "E \t"
       << "index \t"
       << "a' \t"
       << "b' \t"
       << "c' \t"
       << "a \t"
       << "b \t"
       << "c \t" << endl;
  for (int i = 0; i < state.size(); i++) {
    cout << state[i].tz_2 << "\t";
    cout << state[i].n << "\t";
    cout << state[i].l << "\t";
    cout << state[i].j_2 << "\t";
    cout << state[i].E << "\t";
    cout << state[i].index << "\t";
    if (i >= TP_allow[0][0] && i <= TP_allow[1][0]) {
      cout << "*\t";
    } else {
      cout << " \t";
    }
    if (i >= TP_allow[0][1] && i <= TP_allow[1][1]) {
      cout << "*\t";
    } else {
      cout << " \t";
    }
    if (i >= TP_allow[0][2] && i <= TP_allow[1][2]) {
      cout << "*\t";
    } else {
      cout << " \t";
    }
    if (i >= TP_allow[0][3] && i <= TP_allow[1][3]) {
      cout << "*\t";
    } else {
      cout << " \t";
    }
    if (i >= TP_allow[0][4] && i <= TP_allow[1][4]) {
      cout << "*\t";
    } else {
      cout << " \t";
    }
    if (i >= TP_allow[0][5] && i <= TP_allow[1][5]) {
      cout << "*\t";
    } else {
      cout << " \t";
    }
    cout << endl;
  }
}

void SingleP::print_DP(vector<vector<int>> TP_allow) {
  cout << "single particle state" << endl;
  cout << "n_max :" << n_max << "\t l_max" << l_max << endl;
  cout << "tz \t"
       << "n \t"
       << "l \t"
       << "j2 \t"
       << "E \t"
       << "index \t"
       << "a' \t"
       << "b' \t"
       << "a \t"
       << "b \t" << endl;
  for (int i = 0; i < state.size(); i++) {
    cout << state[i].tz_2 << "\t";
    cout << state[i].n << "\t";
    cout << state[i].l << "\t";
    cout << state[i].j_2 << "\t";
    cout << state[i].E << "\t";
    cout << state[i].index << "\t";
    if (i >= TP_allow[0][0] && i <= TP_allow[1][0]) {
      cout << "*\t";
    } else {
      cout << " \t";
    }
    if (i >= TP_allow[0][1] && i <= TP_allow[1][1]) {
      cout << "*\t";
    } else {
      cout << " \t";
    }
    if (i >= TP_allow[0][2] && i <= TP_allow[1][2]) {
      cout << "*\t";
    } else {
      cout << " \t";
    }
    if (i >= TP_allow[0][3] && i <= TP_allow[1][3]) {
      cout << "*\t";
    } else {
      cout << " \t";
    }
    cout << endl;
  }
}

void SingleP::print_DP_Tz(vector<vector<int>> TP_allow) {
  cout << "single particle state" << endl;
  cout << "n_max :" << n_max << "\t l_max" << l_max << endl;
  cout << "tz \t"
       << "n \t"
       << "l \t"
       << "j2 \t"
       << "E \t"
       << "index \t"
       << "a' \t"
       << "b' \t"
       << "a \t"
       << "b \t" << endl;
  for (int i = 0; i < state_Tz.size(); i++) {
    cout << state_Tz[i].tz_2 << "\t";
    cout << state_Tz[i].n << "\t";
    cout << state_Tz[i].l << "\t";
    cout << state_Tz[i].j_2 << "\t";
    cout << state_Tz[i].E << "\t";
    cout << state_Tz[i].index << "\t";
    if (i >= TP_allow[0][0] && i <= TP_allow[1][0]) {
      cout << "*\t";
    } else {
      cout << " \t";
    }
    if (i >= TP_allow[0][1] && i <= TP_allow[1][1]) {
      cout << "*\t";
    } else {
      cout << " \t";
    }
    if (i >= TP_allow[0][2] && i <= TP_allow[1][2]) {
      cout << "*\t";
    } else {
      cout << " \t";
    }
    if (i >= TP_allow[0][3] && i <= TP_allow[1][3]) {
      cout << "*\t";
    } else {
      cout << " \t";
    }

    cout << endl;
  }
}
void SingleP::print(std::string filename) {
  std::ofstream file;
  file.open(filename);
  file << "single particle state" << endl;
  file << "n_max : " << n_max << "\t l_max : " << l_max << endl;
  file << state.size() << " single particle number" << endl;
  file << "tz \t"
       << "n \t"
       << "l \t"
       << "j2 \t"
       << "E \t"
       << "index" << endl;
  for (int i = 0; i < state.size(); i++) {
    file << state[i].tz_2 << "\t";
    file << state[i].n << "\t";
    file << state[i].l << "\t";
    file << state[i].j_2 << "\t";
    file << state[i].E << "\t";
    file << state[i].index << "\t";
    file << endl;
  }
  file.close();
}

void SingleP::print_Tz() {
  cout << "single particle state" << endl;
  cout << "n_max :" << n_max << "\t l_max" << l_max << endl;
  cout << "tz \t"
       << "n \t"
       << "l \t"
       << "j2 \t"
       << "E \t"
       << "index" << endl;
  for (int i = 0; i < state_Tz.size(); i++) {
    cout << state_Tz[i].tz_2 << "\t";
    cout << state_Tz[i].n << "\t";
    cout << state_Tz[i].l << "\t";
    cout << state_Tz[i].j_2 << "\t";
    cout << state_Tz[i].E << "\t";
    cout << state_Tz[i].index << "\t";
    cout << endl;
  }
}
void SingleP::print_Tz(std::string filename) {
  std::ofstream file;
  file.open(filename);
  file << "single particle state" << endl;
  file << "n_max :" << n_max << "\t l_max : " << l_max << endl;
  file << state_Tz.size() << " single particle number" << endl;
  file << "tz \t"
       << "n \t"
       << "l \t"
       << "j2 \t"
       << "E \t"
       << "index" << endl;
  for (int i = 0; i < state_Tz.size(); i++) {
    file << state_Tz[i].tz_2 << "\t";
    file << state_Tz[i].n << "\t";
    file << state_Tz[i].l << "\t";
    file << state_Tz[i].j_2 << "\t";
    file << state_Tz[i].E << "\t";
    file << state_Tz[i].index << "\t";
    file << endl;
  }
  file.close();
}

void SingleP::print_TzM(std::string filename) {
  std::ofstream file;
  file.open(filename);
  file << "single particle state m-scheme" << endl;
  file << "n_max :" << n_max << "\t l_max : " << l_max << endl;
  file << state_TzM.size() << " single particle number" << endl;
  file << "tz \t"
       << "n \t"
       << "l \t"
       << "j2 \t"
       << "m2 \t"
       << "E \t"
       << "index" << endl;
  for (int i = 0; i < state_TzM.size(); i++) {
    file << state_TzM[i].tz_2 << "\t";
    file << state_TzM[i].n << "\t";
    file << state_TzM[i].l << "\t";
    file << state_TzM[i].j_2 << "\t";
    file << state_TzM[i].m_2 << "\t";
    file << state_TzM[i].E << "\t";
    file << state_TzM[i].index << "\t";
    file << endl;
  }
  file.close();
}

int SingleP::FindSP(int n, int l, int j2, int tz_2) {
  int size = state.size();
  int index = -1;
  for (int i = 0; i < size; i++) {
    if (tz_2 == state[i].tz_2 && j2 == state[i].j_2) {
      if (n == state[i].n && l == state[i].l) {
        index = i;
        break;
      }
    }
  }
  return index;
}

//~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~end Single
// Particle~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~//
int DoubleP::Tri(int &j2a, int &j2b, int &j2c) {
  // I+J>=K, I+K>=J, J+K>=I,
  // I/2+J/2+K/2 = INTEGER.
  // TRI=1, WHEN TRIADIC CONDITION IS FULFILLED, TRI=-1 OTHERWISE
  int L2 = j2a + j2b + j2c;
  if (int(L2) != (L2 / 2) * 2)
    return -1;
  L2 = L2 / 2;
  if (j2a * j2b * j2c < 0)
    return -1;
  if ((L2 - j2a) * (L2 - j2b) * (L2 - j2c) < 0)
    return -1;
  return 1;
}

void DoubleP::build_Constrain(SingleP &sp_t) {
  vector<int> abc_vec(4);
  int DP_allow_file_num;
  char u1[256]; // unuse
  std::ifstream DPinput_file;
  DPinput_file.open("data/Input/DP.Input");

  for (int i = 0; i < 4; i++) {

    int A = -1;
    DPinput_file >> A;
    // cout << A << "\t** " << endl;
    abc_vec[i] = A;
    DPinput_file.getline(u1, 256, '\n');
  }
  DPinput_file >> DP_allow_file_num;
  DPinput_file.close();

  int a_min = abc_vec[0];
  int a_max = abc_vec[1];
  int b_min = abc_vec[2];
  int b_max = abc_vec[3];

  vector<vector<vector<int>>> DP_allow_gather_Tz;
  vector<int> E_cut_gather;
  vector<int> Same_Particles_gather;
  vector<vector<int>> E_sign_gather;

  //
  // int TP_allow_file_num = 2;

  std::ifstream DP_Allow_file;
  for (int i = 0; i < DP_allow_file_num; i++) {
    vector<vector<int>> DP_allow_vec;
    vector<vector<int>> DP_allow_vec_Tz;

    DP_allow_vec.resize(2);
    DP_allow_vec_Tz.resize(2);
    for (int j = 0; j < 2; j++) {
      DP_allow_vec[j].resize(4);
      DP_allow_vec_Tz[j].resize(4);
    }
    vector<int> E_sign(4);
    int Same_Particles_Flag;
    int E_cut;
    int ii = i;
    std::string filename = "data/Input/DP_allow_" + std::to_string(ii);
    filename += ".Input";
    DP_Allow_file.open(filename);
    DP_Allow_file >> E_cut;
    DP_Allow_file.getline(u1, 256, '\n');
    DP_Allow_file >> Same_Particles_Flag;
    DP_Allow_file.getline(u1, 256, '\n');
    for (int i = 0; i < 4; i++) {
      int min, max;
      DP_Allow_file >> min >> max >> E_sign[i];
      // cout << "\t TP allow : " << min << "\t" << max << endl;
      // cout << A << "\t** " << endl;
      DP_allow_vec[0][i] = min;
      DP_allow_vec[1][i] = max;
      DP_allow_vec_Tz[0][i] = min * 2;
      DP_allow_vec_Tz[1][i] = max * 2 + 1;
      DP_Allow_file.getline(u1, 256, '\n');
    }
    DP_Allow_file.close();
    DP_allow_gather_Tz.push_back(DP_allow_vec_Tz);
    E_cut_gather.push_back(E_cut);
    E_sign_gather.push_back(E_sign);
    Same_Particles_gather.push_back(Same_Particles_Flag);

    sp_t.print_DP(DP_allow_vec);
  }
  // cout<<"a_max : "<<a_max<<"\t b_max : "<<b_max<<endl;
  this->build_Tz(sp_t, 2 * a_min, 2 * a_max + 1, 2 * b_min, 2 * b_max + 1);
  // this->print_Tz("data/DP_Tz_un1.dat");
  // this->print_channel_Tz("data/DP_channel_Tz_un1.dat");
  this->Select_FI_Tz(sp_t, DP_allow_gather_Tz);
  // this->print_Tz("data/DP_Tz_un2.dat");
  // this->print_channel_Tz("data/DP_channel_Tz_un2.dat");
  this->Delet_Sort_Tz();

  this->print_Tz("data/DP_Tz.dat");
  this->print_channel_Tz("data/DP_channel_Tz.dat");
}

void DoubleP::build_Tz(SingleP &sp_t, int a_beg, int a_end, int b_beg,
                       int b_end) {
  int j2a_max_a = -1;
  int j2a_max_b = -1;

  int size_sp = sp_t.state_Tz.size();
  if (a_end > size_sp || b_end > size_sp) {
    cout << "Wrong @ DoubleP::build_Tz " << endl;
    exit(0);
  }

  for (int a = a_beg; a <= a_end; a++) {
    int j2_t = sp_t.state_Tz[a].j_2;
    if (j2_t > j2a_max_a)
      j2a_max_a = j2_t;
  }
  for (int b = b_beg; b <= b_end; b++) {
    int j2_t = sp_t.state_Tz[b].j_2;
    if (j2_t > j2a_max_b)
      j2a_max_b = j2_t;
  }

  J2_max = j2a_max_a + j2a_max_b;
  int sp_size = sp_t.state_Tz.size();
  int chan_size = (J2_max) / 2 + 1;
  chan_size = chan_size * 3; // Tz = -1,0,1
  state_Tz.clear();
  state_Tz.resize(chan_size);
  cout << "DP size_channel : " << chan_size << endl;
  /*cout << "Test size = " << sp_size << "\t ch_size = " << chan_size << endl;
  cout << a_beg << "\t" << a_end << "\t\t" << b_beg << "\t" << b_end << "\t\t"
       << c_end << "\t" << c_end << endl;/**/
  State_DP state_t;
  int tot_num = 0;
  int chan_num = -1;
  for (int T2z = -2; T2z <= 2; T2z = T2z + 2) {
    for (int J_2 = 0; J_2 <= J2_max; J_2 = J_2 + 2) {
      chan_num++;
      // cout<<"chan : "<<chan_num<<"\t T2z : "<<T2z<<"\t J_2 : "<<J_2<<endl;
      for (int a = a_beg; a <= a_end; a++) {
        for (int b = b_beg; b <= b_end; b++) {
          if (T2z != (sp_t.state_Tz[a].tz_2 + sp_t.state_Tz[b].tz_2))
            continue;
          if (Tri(J_2, sp_t.state_Tz[a].j_2, sp_t.state_Tz[b].j_2) != 1)
            continue;
          // if (sp_t.state_Tz[a].l == 0 && sp_t.state_Tz[b].l == 0 && a == b)
          //   continue;
          state_t.I = -1;
          state_t.F = -1;
          state_t.a = a;
          state_t.b = b;
          state_t.J_2 = J_2;
          state_t.T_2z = T2z;
          state_t.E = sp_t.state_Tz[a].E + sp_t.state_Tz[b].E;
          this->state_Tz[chan_num].push_back(state_t);
          tot_num++;
        }
      }
    }
  }
  cout << "DP total Num : " << tot_num << endl;
  this->sort_Tz();
}

void DoubleP::sort_Tz() {

  // Sort state
  auto it = state_Tz.begin();
  while (it != state_Tz.end()) {
    if (it->size() == 0) {
      it = state_Tz.erase(it);
    } else {
      it++;
    }
  }
  int chan_size = state_Tz.size();
  for (int i = 0; i < chan_size; i++) {
    stable_sort(state_Tz[i].begin(), state_Tz[i].end(),
                [](State_DP a, State_DP b) {
                  return (1000000 * a.E + 100000 * a.T_2z + 10000 * a.J_2 +
                          100 * a.a + a.b) < (1000000 * b.E + 100000 * b.T_2z +
                                              10000 * b.J_2 + 100 * b.a + b.b);
                });
    int index = 0;
    for (int j = 0; j < state_Tz[i].size(); j++) {
      state_Tz[i][j].index = index;
      index++;
    }
  }

  // cout<<"good? 1 "<<endl;
  // Rebuild Channel
  chan_size = state_Tz.size();
  channel_Tz.clear();

  int all_size_T_t = 0;
  for (int i = 0; i < chan_size; i++) {
    all_size_T_t += state_Tz[i].size();
  }
  this->all_size_Tz = all_size_T_t;

  channel_Tz.resize(chan_size);

  for (int i = 0; i < chan_size; i++) {
    if (state_Tz[i].size() == 0) {
      channel_Tz[i].num = 0;
      channel_Tz[i].J_2 = 0;
      channel_Tz[i].T_2 = 0;
      channel_Tz[i].E_min = 0;
      channel_Tz[i].E_max = 0;
      // channel.erase(i);
    } else {
      channel_Tz[i].num = state_Tz[i].size();
      channel_Tz[i].J_2 = state_Tz[i][0].J_2;
      channel_Tz[i].T_2z = state_Tz[i][0].T_2z;
      int E_min = state_Tz[i][0].E;
      int E_max = state_Tz[i][0].E;
      for (int j = 0; j < state_Tz[i].size(); j++) {
        if (E_min > state_Tz[i][j].E)
          E_min = state_Tz[i][j].E;
        if (E_max < state_Tz[i][j].E)
          E_max = state_Tz[i][j].E;
      }
      channel_Tz[i].E_min = E_min;
      channel_Tz[i].E_max = E_max;
    }
  }
  // cout<<"good? 2 "<<endl;

  // mark the indexE
  for (int i = 0; i < chan_size; i++) {
    int Emin = channel_Tz[i].E_min;
    int Emax = channel_Tz[i].E_max;
    channel_Tz[i].indexE.resize(Emax + 1);
    channel_Tz[i].indexE_beg.resize(Emax + 1);
    int flag = 1;
    for (int j = 0; j < state_Tz[i].size(); j++) {
      channel_Tz[i].indexE[state_Tz[i][j].E] = state_Tz[i][j].index;
    }
    for (int E = Emin; E <= Emax; E++) {
      int indexE_beg = 0;
      for (int j = 0; j < state_Tz[i].size(); j++) {
        if (state_Tz[i][j].E == E) {
          indexE_beg = j;
          break;
        }
      }
      channel_Tz[i].indexE_beg[E] = indexE_beg;
    }
  }
  // mark channel_TP.E_min/max_F/I
  for (int ch = 0; ch < chan_size; ch++) {
    int Emin_f_t = 100;
    int Emin_i_t = 100;
    int Emax_f_t = 0;
    int Emax_i_t = 0;
    int Emax_t = 0;
    int Emin_t = 100;

    for (int i = 0; i < state_Tz[ch].size(); i++) {
      int E_t = state_Tz[ch][i].E;
      if (Emin_t > E_t)
        Emin_t = E_t;
      if (Emax_t < E_t)
        Emax_t = E_t;
      if (state_Tz[ch][i].F > 0) {
        Emin_f_t = Emin_t;
        Emax_f_t = Emax_t;
      }
      if (state_Tz[ch][i].I > 0) {
        Emin_i_t = Emin_t;
        Emax_i_t = Emax_t;
      }
    }

    this->channel_Tz[ch].E_min_F = Emin_f_t;
    this->channel_Tz[ch].E_min_I = Emin_i_t;
    this->channel_Tz[ch].E_max_F = Emax_f_t;
    this->channel_Tz[ch].E_max_I = Emax_i_t;
  }
  // cout<<"good? 3 "<<endl;

  // Find J2_max
  int J2_max_t = 1;
  for (int i = 0; i < chan_size; i++) {
    if (channel_Tz[i].J_2 > J2_max_t)
      J2_max_t = channel_Tz[i].J_2;
  }
  // cout<<"good? 4 "<<endl;

  this->J2_max_DP = J2_max_t;

  // cout<<"good? 5 "<<endl;
  // cout<<"good? 6 "<<endl;
}
int DoubleP::FindC_Tz(int J2, int T_2z) {
  // Find Channel with J2 T2
  int res = -1;
  int size = channel_Tz.size();
  if (size == 0)
    return -1;
  for (int i = 0; i < size; i++) {
    if (J2 == channel_Tz[i].J_2 && T_2z == channel_Tz[i].T_2z) {
      res = i;
      break;
    }
  }
  if (res == -1) {
    cout << "wrong at find channel in DoubleP::FindC_Tz " << endl;
    cout << "\t J2 = " << J2 << "\t T_2z = " << T_2z << endl;
    exit(0);
  }
  return res;
}

int DoubleP::FindDP_Tz(int J2, int T_2z, int a, int b) {
  int index = -1;
  int ch_Tz = FindC_Tz(J2, T_2z);
  int size = state_Tz[ch_Tz].size();
  for (int i = 0; i < size; i++) {
    int a_t = state_Tz[ch_Tz][i].a;
    int b_t = state_Tz[ch_Tz][i].b;

    if (b_t != b || a_t != a)
      continue;
    index = i;
    break;
  }
  if (index == -1) {
    cout << "wrong at find DP in DoubleP::FindDP_Tz " << endl;
    cout << "J2 = " << J2 << "\t T_2z = " << T_2z << "\t ch_Tz = " << ch_Tz
         << "\t a = " << a << "\t b = " << b << endl;
    exit(0);
  }
  return index;
}

void DoubleP::print_Tz(std::string filename) {
  std::ofstream file;
  file.open(filename);
  file << "a b JJ-coupling state_Tz" << endl;
  file << "Energy tot Max is : " << Emax << " ( hbar-omega )" << endl;
  file << "total number of basis is : " << endl;

  file << all_size_Tz << endl;
  file << channel_Tz.size() << "\t channel_Tz number" << endl;
  for (int i = 0; i < channel_Tz.size(); i++) {
    file << state_Tz[i].size() << "\t";
  }
  file << endl;
  for (int i = 0; i < channel_Tz.size(); i++) {
    file << "===== J channel_Tz : =====" << i << endl;
    file << "J2 \t"
         << "T_2z \t"
         << "a \t"
         << "b \t"
         << "E \t"
         << "index \t F\t I" << endl;
    for (int j = 0; j < state_Tz[i].size(); j++) {
      file << state_Tz[i][j].J_2 << "\t";
      file << state_Tz[i][j].T_2z << "\t";
      file << state_Tz[i][j].a << "\t";
      file << state_Tz[i][j].b << "\t";
      file << state_Tz[i][j].E << "\t";
      file << state_Tz[i][j].index << "\t";
      file << state_Tz[i][j].F << "\t";
      file << state_Tz[i][j].I << "\t";
      /*if (state_Tz[i][j].F > 0)
        file << "*";
      file << "\t";
      if (state_Tz[i][j].I > 0)
        file << "* \t";/**/

      file << endl;
    }
  }
  file.close();
}

void DoubleP::print_channel_Tz(std::string filename) {
  std::ofstream file;
  file.open(filename);
  int size = this->channel_Tz.size();
  file << " DP channel_Tz Total Num: " << size << endl;
  file << "index \t"
       << "num \t"
       << "J2 \t"
       << "T_2z \t"
       << "E_min \t"
       << "E_max \t" << endl;
  for (int i = 0; i < size; i++) {
    file << i << "\t" << channel_Tz[i].num << "\t" << channel_Tz[i].J_2 << "\t"
         << channel_Tz[i].T_2z << "\t" << channel_Tz[i].E_min << "\t"
         << channel_Tz[i].E_max << endl;
  }
  file << "E_min and E_max index:" << endl;
  file << "\t \t   "
       << "E : 0\t1\t2 \t3 \t4 \t5 ....." << endl;
  for (int i = 0; i < size; i++) {
    file << "channel : " << i << endl;
    file << "\t E_min index : ";
    for (int j = 0; j < channel_Tz[i].indexE_beg.size(); j++) {
      file << channel_Tz[i].indexE_beg[j] << "\t";
    }
    file << endl;
    file << "\t E_max index : ";
    for (int j = 0; j < channel_Tz[i].indexE.size(); j++) {
      file << channel_Tz[i].indexE[j] << "\t";
    }
    file << endl;
  }
  file.close();
}

void DoubleP::Select_FI_Tz(SingleP &sp_t,
                           vector<vector<vector<int>>> allow_FI_gather) {
  // orbits which particles can run :
  //  min: allow_FI[0][a', b', c', a, b, c]
  //  max: allow_FI[1][a', b', c', a, b, c]
  int file_num = allow_FI_gather.size();
  int sp_min_t = 100;
  int sp_max_t = -1;

  int Flag = 0;
  for (int k = 0; k < file_num; k++) {
    Flag = std::pow(10, k);

    vector<vector<int>> allow_FI;
    allow_FI = allow_FI_gather[k];
    // cout << "\ta'\t b'\t c'\t a\t b\t c\t" << endl;
    /*for (int i = 0; i < 2; i++) {
      if (i == 0)
        cout << "min: \t";
      if (i == 1)
        cout << "max: \t";
      for (int j = 0; j < 4; j++) {
        cout << allow_FI[i][j] << "\t";
        if (allow_FI[i][j] > sp_max_t)
          sp_max_t = allow_FI[i][j];
        if (allow_FI[i][j] < sp_min_t)
          sp_min_t = allow_FI[i][j];
      }
      cout << endl;
    }/**/

    int ap_min = allow_FI[0][0];
    int bp_min = allow_FI[0][1];
    // int cp_min = allow_FI[0][2];
    int a_min = allow_FI[0][2];
    int b_min = allow_FI[0][3];
    // int c_min = allow_FI[0][5];

    int ap_max = allow_FI[1][0];
    int bp_max = allow_FI[1][1];
    // int cp_max = allow_FI[1][2];
    int a_max = allow_FI[1][2];
    int b_max = allow_FI[1][3];
    // int c_max = allow_FI[1][5];

    for (int ch = 0; ch < state_Tz.size(); ch++) {
      for (int i = 0; i < state_Tz[ch].size(); i++) {
        int a = state_Tz[ch][i].a;
        int b = state_Tz[ch][i].b;
        // int c = state[ch][i].c;
        if (a < ap_min || a > ap_max || b < bp_min || b > bp_max) {
          int unuse;
        } else {
          if (state_Tz[ch][i].F < 0) {
            state_Tz[ch][i].F = Flag;
          } else {
            state_Tz[ch][i].F += Flag;
          }
        }
        if (a < a_min || a > a_max || b < b_min || b > b_max) {
          int unuse;
        } else {
          if (state_Tz[ch][i].I < 0) {
            state_Tz[ch][i].I = Flag;
          } else {
            state_Tz[ch][i].I += Flag;
          }
        }
      }
    }
  }
  // Find j2max_sp
  int j2max_sp_t = 1;
  int sp_min = sp_min_t;
  int sp_max = sp_max_t;

  for (int i = sp_min; i <= sp_max; i++) {
    int j2 = sp_t.state_Tz[i].j_2;
    if (j2 > j2max_sp_t)
      j2max_sp_t = j2;
  }
  this->j2max_sp = j2max_sp_t;
}

void DoubleP::Delet_Sort_Tz() {
  // call Select_FI_T first !!!
  // 1. Delete by .F = -1 && .I = -1
  // 2. Sort for index and channel
  //
  // Delete
  vector<vector<State_DP>> state_Tz_t;
  state_Tz_t = this->state_Tz;
  int chan_size = state_Tz_t.size();
  state_Tz.clear();
  state_Tz.resize(chan_size);

  for (int ch = 0; ch < chan_size; ch++) {
    int num_I = 0;
    int num_F = 0;
    for (int i = 0; i < state_Tz_t[ch].size(); i++) {
      if (state_Tz_t[ch][i].F > 0)
        num_F++;
      if (state_Tz_t[ch][i].I > 0)
        num_I++;
    }
    if (num_F > 0 && num_I > 0) {
      for (int i = 0; i < state_Tz_t[ch].size(); i++) {
        if (state_Tz_t[ch][i].F < 0 && state_Tz_t[ch][i].I < 0)
          continue;
        state_Tz[ch].push_back(state_Tz_t[ch][i]);
      }
    }
  }
  this->sort_Tz();
}
//~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~start Triple
// Particle~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~//
void TripleP::build_Tz(SingleP &sp_t) {
  int j2a_max = sp_t.l_max * 2 + 1;
  J2_max = 3 * j2a_max;
  int sp_size = sp_t.state_Tz.size();
  int chan_size = J2_max + 1;
  state_Tz.resize(chan_size);

  State_TP state_t;

  int chan_num = -1;
  for (int T_2 = 1; T_2 <= 3; T_2 = T_2 + 2) {
    for (int J_2 = 1; J_2 <= J2_max; J_2 = J_2 + 2) {
      chan_num++;
      int T2ab_min = 2;
      int T2ab_max = 2;
      if (T_2 == 1) {
        T2ab_min = 0;
      }
      int J2ab_min, J2ab_max;
      J2ab_min = std::max(0, J_2 - j2a_max);
      J2ab_max = J_2 + j2a_max;
      for (int T2ab = T2ab_min; T2ab <= T2ab_max; T2ab++) {
        for (int J2ab = J2ab_min; J2ab <= J2ab_max; J2ab = J2ab + 2) {
          for (int a = 0; a < sp_size; a++) {
            for (int b = 0; b < sp_size; b++) {
              // for(int b=0;b<1;b++){
              if (T2ab !=
                  std::abs(sp_t.state_Tz[a].tz_2 + sp_t.state_Tz[b].tz_2))
                continue;
              if (Tri(J2ab, sp_t.state_Tz[a].j_2, sp_t.state_Tz[b].j_2) != 1)
                continue;
              for (int c = 0; c < sp_size; c++) {
                // for(int c=0;c<1;c++){
                if (std::abs(T2ab + sp_t.state_Tz[c].tz_2) != T_2)
                  continue;
                if (Tri(J2ab, sp_t.state_Tz[c].j_2, J_2) != 1)
                  continue;
                if (sp_t.state_Tz[a].l == 0 && sp_t.state_Tz[b].l == 0 &&
                    a == b && b == c)
                  continue;
                state_t.a = a;
                state_t.b = b;
                state_t.c = c;
                state_t.Jab_2 = J2ab;
                state_t.Tab_2z = T2ab;
                state_t.J_2 = J_2;
                state_t.T_2z = T_2;
                state_t.E = sp_t.state_Tz[a].E + sp_t.state_Tz[b].E +
                            sp_t.state_Tz[c].E;
                state_t.Par = this->phase(state_t.E);
                state_Tz[chan_num].push_back(state_t);
              }
            }
          }
        }
      }
    }
  }

  for (int i = 0; i < chan_size; i++) {
    stable_sort(state_Tz[i].begin(), state_Tz[i].end(),
                [](State_TP a, State_TP b) {
                  return (100000 * a.Tab_2 + 10000 * a.Jab_2 + 1000 * a.E +
                          100 * a.a + a.b + 0.01 * a.c) <
                         (100000 * b.Tab_2 + 10000 * b.Jab_2 + 1000 * b.E +
                          100 * b.a + b.b + 0.01 * b.c);
                });
    int index = 0;
    for (int j = 0; j < state_Tz[i].size(); j++) {
      state_Tz[i][j].index = index;
      index++;
    }
  }

  channel_Tz.resize(chan_size);

  for (int i = 0; i < chan_size; i++) {
    if (state_Tz[i].size() == 0) {
      channel_Tz[i].num = 0;
      channel_Tz[i].J_2 = 0;
      channel_Tz[i].T_2 = 0;
      channel_Tz[i].E_min = 0;
      channel_Tz[i].E_max = 0;
    } else {
      channel_Tz[i].num = state_Tz[i].size();
      channel_Tz[i].J_2 = state_Tz[i][0].J_2;
      channel_Tz[i].T_2z = state_Tz[i][0].T_2z;
      int E_min = state_Tz[i][0].E;
      int E_max = state_Tz[i][0].E;
      for (int j = 0; j < state_Tz[i].size(); j++) {
        if (E_min > state_Tz[i][j].E)
          E_min = state_Tz[i][j].E;
        if (E_max < state_Tz[i][j].E)
          E_max = state_Tz[i][j].E;
      }
      channel_Tz[i].E_min = E_min;
      channel_Tz[i].E_max = E_max;
    }
  }

  // calulate JT2channel_Tz
}

void TripleP::build_Tz(SingleP &sp_t, int Emax_t) {
  Emax = Emax_t;
  int j2a_max = sp_t.l_max * 2 + 1;
  J2_max = 3 * j2a_max;
  int sp_size = sp_t.state_Tz.size();
  int chan_size = (J2_max + 1) / 2;
  state_Tz.resize(chan_size);

  State_TP state_t;

  int chan_num = -1;
  for (int J_2 = 1; J_2 <= J2_max; J_2 = J_2 + 2) {
    chan_num++;
    for (int T2z = -3; T2z <= 3; T2z = T2z + 2) {
      int T2abz_min = -2;
      int T2abz_max = 2;

      int J2ab_min, J2ab_max;
      J2ab_min = std::max(0, J_2 - j2a_max);
      J2ab_max = J_2 + j2a_max;
      for (int T2abz = T2abz_min; T2abz <= T2abz_max; T2abz = T2abz + 2) {
        for (int J2ab = J2ab_min; J2ab <= J2ab_max; J2ab = J2ab + 2) {
          for (int a = 0; a < sp_size; a++) {
            if (sp_t.state_Tz[a].E > Emax)
              continue;
            for (int b = 0; b < sp_size; b++) {
              // for(int b=0; b<1; b++){
              if (sp_t.state_Tz[a].E + sp_t.state_Tz[b].E > Emax)
                continue;
              if (T2abz != (sp_t.state_Tz[a].tz_2 + sp_t.state_Tz[b].tz_2))
                continue;
              if (Tri(J2ab, sp_t.state_Tz[a].j_2, sp_t.state_Tz[b].j_2) != 1)
                continue;
              for (int c = 0; c < sp_size; c++) {
                // for(int c=0; c<1; c++){
                if (sp_t.state_Tz[a].E + sp_t.state_Tz[b].E +
                        sp_t.state_Tz[c].E >
                    Emax)
                  continue;
                if ((sp_t.state_Tz[a].tz_2 + sp_t.state_Tz[b].tz_2 +
                     sp_t.state_Tz[c].tz_2) != T2z)
                  continue;
                if (Tri(J2ab, sp_t.state_Tz[c].j_2, J_2) != 1)
                  continue;
                if (sp_t.state_Tz[a].l == 0 && sp_t.state_Tz[b].l == 0 &&
                    a == b && b == c)
                  continue;
                state_t.a = a;
                state_t.b = b;
                state_t.c = c;
                state_t.Jab_2 = J2ab;
                state_t.Tab_2z = T2abz;
                state_t.J_2 = J_2;
                state_t.T_2z = T2z;
                state_t.E = sp_t.state_Tz[a].E + sp_t.state_Tz[b].E +
                            sp_t.state_Tz[c].E;
                state_t.Par = this->phase(state_t.E);
                state_Tz[chan_num].push_back(state_t);
              }
            }
          }
        }
      }
    }
  }

  for (int i = 0; i < chan_size; i++) {
    stable_sort(state_Tz[i].begin(), state_Tz[i].end(),
                [](State_TP a, State_TP b) {
                  return (100000 * a.E + 10000 * a.Jab_2 + 100 * a.a + a.b +
                          0.01 * a.c + 0.001 * a.Tab_2z + 0.0001 * a.T_2z) <
                         (100000 * b.E + 10000 * b.Jab_2 + 100 * b.a + b.b +
                          0.01 * b.c + 0.001 * b.Tab_2z + 0.0001 * b.T_2z);
                });
    int index = 0;
    for (int j = 0; j < state_Tz[i].size(); j++) {
      state_Tz[i][j].index = index;
      index++;
    }
  }

  channel_Tz.resize(chan_size);

  for (int i = 0; i < chan_size; i++) {
    if (state_Tz[i].size() == 0) {
      channel_Tz[i].num = 0;
      channel_Tz[i].J_2 = 0;
      channel_Tz[i].T_2z = 0;
      channel_Tz[i].E_min = 0;
      channel_Tz[i].E_max = 0;
    } else {
      channel_Tz[i].num = state_Tz[i].size();
      channel_Tz[i].J_2 = state_Tz[i][0].J_2;
      channel_Tz[i].T_2z = state_Tz[i][0].T_2z;
      int E_min = state_Tz[i][0].E;
      int E_max = state_Tz[i][0].E;
      for (int j = 0; j < state_Tz[i].size(); j++) {
        if (E_min > state_Tz[i][j].E)
          E_min = state_Tz[i][j].E;
        if (E_max < state_Tz[i][j].E)
          E_max = state_Tz[i][j].E;
      }
      channel_Tz[i].E_min = E_min;
      channel_Tz[i].E_max = E_max;
    }
  }

  // mark the indexE
  for (int i = 0; i < chan_size; i++) {
    int Emin = channel_Tz[i].E_min;
    int Emax = channel_Tz[i].E_max;
    channel_Tz[i].indexE.resize(Emax + 1);
    channel_Tz[i].indexE_beg.resize(Emax + 1);
    int flag = 1;
    for (int j = 0; j < state_Tz[i].size(); j++) {
      channel_Tz[i].indexE[state_Tz[i][j].E] = state_Tz[i][j].index;
    }
    for (int E = Emin; E <= Emax; E++) {
      int indexE_beg = 0;
      for (int j = 0; j < state_Tz[i].size(); j++) {
        if (state_Tz[i][j].E == E) {
          indexE_beg = j;
          break;
        }
      }
      channel_Tz[i].indexE_beg[E] = indexE_beg;
    }
  }
}

void TripleP::build_Tz(SingleP &sp_t, int a_beg, int a_end, int b_beg,
                       int b_end, int c_beg, int c_end) {

  // Emax = Emax_t;
  int j2a_max_a = 1;
  int j2a_max_b = 1;
  int j2a_max_c = 1;

  for (int a = a_beg; a <= a_end; a++) {
    int j2_t = sp_t.state_Tz[a].j_2;
    if (j2_t > j2a_max_a)
      j2a_max_a = j2_t;
  }
  for (int b = b_beg; b <= b_end; b++) {
    int j2_t = sp_t.state_Tz[b].j_2;
    if (j2_t > j2a_max_b)
      j2a_max_b = j2_t;
  }
  for (int c = c_beg; c <= c_end; c++) {
    int j2_t = sp_t.state_Tz[c].j_2;
    if (j2_t > j2a_max_c)
      j2a_max_c = j2_t;
  }

  J2_max = j2a_max_a + j2a_max_b + j2a_max_c;
  int sp_size = sp_t.state_Tz.size();
  int chan_size = (J2_max + 1) / 2;
  chan_size = chan_size * 4; // Tz = -3,-1,1,3
  state_Tz.resize(chan_size);
  /*cout << "Test size = " << sp_size << "\t ch_size = " << chan_size << endl;
  cout << a_beg << "\t" << a_end << "\t\t" << b_beg << "\t" << b_end << "\t\t"
       << c_end << "\t" << c_end << endl;/**/
  State_TP state_t;

  int chan_num = -1;
  for (int T2z = -3; T2z <= 3; T2z = T2z + 2) {
    for (int J_2 = 1; J_2 <= J2_max; J_2 = J_2 + 2) {
      chan_num++;

      int T2abz_min = -2;
      int T2abz_max = 2;

      int J2ab_min, J2ab_max;
      J2ab_min = 0;
      J2ab_max = J_2 + j2a_max_c;
      for (int T2abz = T2abz_min; T2abz <= T2abz_max; T2abz = T2abz + 2) {
        for (int J2ab = J2ab_min; J2ab <= J2ab_max; J2ab = J2ab + 2) {
          for (int a = a_beg; a <= a_end; a++) {
            for (int b = b_beg; b <= b_end; b++) {
              if (T2abz != (sp_t.state_Tz[a].tz_2 + sp_t.state_Tz[b].tz_2))
                continue;
              if (Tri(J2ab, sp_t.state_Tz[a].j_2, sp_t.state_Tz[b].j_2) != 1)
                continue;
              for (int c = c_beg; c <= c_end; c++) {
                if ((sp_t.state_Tz[a].tz_2 + sp_t.state_Tz[b].tz_2 +
                     sp_t.state_Tz[c].tz_2) != T2z)
                  continue;
                if (Tri(J2ab, sp_t.state_Tz[c].j_2, J_2) != 1)
                  continue;
                if (sp_t.state_Tz[a].l == 0 && sp_t.state_Tz[b].l == 0 &&
                    a == b && b == c)
                  continue;

                state_t.I = -1;
                state_t.F = -1;
                state_t.Ecut_Flag = -1;
                state_t.a = a;
                state_t.b = b;
                state_t.c = c;
                state_t.Jab_2 = J2ab;
                state_t.Tab_2z = T2abz;
                state_t.J_2 = J_2;
                state_t.T_2z = T2z;
                state_t.E = sp_t.state_Tz[a].E + sp_t.state_Tz[b].E +
                            sp_t.state_Tz[c].E;
                state_t.Par = this->phase(state_t.E);
                state_Tz[chan_num].push_back(state_t);
              }
            }
          }
        }
      }
    }
  }

  for (int i = 0; i < chan_size; i++) {
    stable_sort(state_Tz[i].begin(), state_Tz[i].end(),
                [](State_TP a, State_TP b) {
                  return (100000 * a.E + 10000 * a.Jab_2 + 100 * a.a + a.b +
                          0.01 * a.c + 0.001 * a.Tab_2z + 0.0001 * a.T_2z) <
                         (100000 * b.E + 10000 * b.Jab_2 + 100 * b.a + b.b +
                          0.01 * b.c + 0.001 * b.Tab_2z + 0.0001 * b.T_2z);
                });
    int index = 0;
    for (int j = 0; j < state_Tz[i].size(); j++) {
      state_Tz[i][j].index = index;
      index++;
    }
  }

  // aLL_size
  int all_size_Tz_t = 0;
  for (int i = 0; i < chan_size; i++) {
    all_size_Tz += chan_size;
  }
  this->all_size_Tz = all_size_Tz_t;

  channel_Tz.resize(chan_size);

  for (int i = 0; i < chan_size; i++) {
    if (state_Tz[i].size() == 0) {
      channel_Tz[i].num = 0;
      channel_Tz[i].J_2 = 0;
      channel_Tz[i].T_2z = 0;
      channel_Tz[i].E_min = 0;
      channel_Tz[i].E_max = 0;
    } else {
      channel_Tz[i].num = state_Tz[i].size();
      channel_Tz[i].J_2 = state_Tz[i][0].J_2;
      channel_Tz[i].T_2z = state_Tz[i][0].T_2z;
      int E_min = state_Tz[i][0].E;
      int E_max = state_Tz[i][0].E;
      for (int j = 0; j < state_Tz[i].size(); j++) {
        if (E_min > state_Tz[i][j].E)
          E_min = state_Tz[i][j].E;
        if (E_max < state_Tz[i][j].E)
          E_max = state_Tz[i][j].E;
      }
      channel_Tz[i].E_min = E_min;
      channel_Tz[i].E_max = E_max;
    }
  }

  // mark the indexE
  for (int i = 0; i < chan_size; i++) {
    int Emin = channel_Tz[i].E_min;
    int Emax = channel_Tz[i].E_max;
    channel_Tz[i].indexE.resize(Emax + 1);
    channel_Tz[i].indexE_beg.resize(Emax + 1);
    int flag = 1;
    for (int j = 0; j < state_Tz[i].size(); j++) {
      channel_Tz[i].indexE[state_Tz[i][j].E] = state_Tz[i][j].index;
    }
    for (int E = Emin; E <= Emax; E++) {
      int indexE_beg = 0;
      for (int j = 0; j < state_Tz[i].size(); j++) {
        if (state_Tz[i][j].E == E) {
          indexE_beg = j;
          break;
        }
      }
      channel_Tz[i].indexE_beg[E] = indexE_beg;
    }
  } /**/
}

void TripleP::build_TzM(SingleP &SP, int Emax_t) {
  Emax = Emax_t;
  State_TPM state_t;
  int size = SP.state_TzM.size();
  for (int a = 0; a < size; a++) {
    for (int b = 0; b < size; b++) {
      if (b == a)
        continue;
      int Tab_2z = SP.state_TzM[a].tz_2 + SP.state_TzM[b].tz_2;
      int Jab_2m = SP.state_TzM[a].m_2 + SP.state_TzM[b].m_2;
      for (int c = 0; c < size; c++) {
        if (c == b || c == a)
          continue;
        int T_2z = Tab_2z + SP.state_TzM[c].tz_2;
        int J_2m = Jab_2m + SP.state_TzM[c].m_2;
        int E = SP.state_TzM[a].E + SP.state_TzM[b].E + SP.state_TzM[c].E;
        state_t.a_m = a;
        state_t.b_m = b;
        state_t.c_m = c;
        state_t.Jab_2m = Jab_2m;
        state_t.Tab_2z = Tab_2z;
        state_t.J_2m = J_2m;
        state_t.T_2z = T_2z;
        state_t.E = E;
        state_TzM.push_back(state_t);
      }
    }
  }
  stable_sort(state_TzM.begin(), state_TzM.end(), [](State_TPM a, State_TPM b) {
    return (1000000 * a.E + 100000 * a.T_2z + 10000 * a.J_2m + 1000 * a.a_m +
            10 * a.b_m + 0.1 * a.c_m + 0.01 * a.Tab_2z + 0.001 * a.Jab_2m) <
           (1000000 * b.E + 100000 * b.T_2z + 10000 * b.J_2m + 1000 * b.a_m +
            10 * b.b_m + 0.1 * b.c_m + 0.01 * b.Tab_2z + 0.001 * b.Jab_2m);
  });
  for (int i = 0; i < state_TzM.size(); i++) {
    state_TzM[i].index = i;
  }
}

void TripleP::build_T(SingleP &sp_t, int Emax_t) {
  Emax = Emax_t;
  int j2a_max = sp_t.l_max * 2 + 1;
  J2_max = 3 * j2a_max;
  int sp_size = sp_t.state.size();
  int chan_size = J2_max + 1;
  state.resize(chan_size);

  State_TP state_t;

  int chan_num = -1;
  for (int T_2 = 1; T_2 <= 3; T_2 = T_2 + 2) {
    for (int J_2 = 1; J_2 <= J2_max; J_2 = J_2 + 2) {
      chan_num++;
      int T2ab_min = 2;
      int T2ab_max = 2;
      if (T_2 == 1) {
        T2ab_min = 0;
      }
      int J2ab_min, J2ab_max;
      J2ab_min = std::max(0, J_2 - j2a_max);
      J2ab_max = J_2 + j2a_max;
      for (int T2ab = T2ab_min; T2ab <= T2ab_max; T2ab = T2ab + 2) {
        for (int T_2z = -T_2; T_2z <= T_2; T_2z = T_2z + 2) {
          // cout<<T_2z<<endl;
          for (int J2ab = J2ab_min; J2ab <= J2ab_max; J2ab = J2ab + 2) {
            for (int a = 0; a < sp_size; a++) {
              // for(int a=0; a<1; a++){
              if (sp_t.state[a].E > Emax)
                continue;
              for (int b = 0; b < sp_size; b++) {
                // for(int b=0; b<1; b++){
                if (sp_t.state[a].E + sp_t.state[b].E > Emax)
                  continue;
                if (Tri(J2ab, sp_t.state[a].j_2, sp_t.state[b].j_2) != 1)
                  continue;
                for (int c = 0; c < sp_size; c++) {
                  // for(int c=0; c<1; c++){
                  if (sp_t.state[a].E + sp_t.state[b].E + sp_t.state[c].E >
                      Emax)
                    continue;

                  if (Tri(J2ab, sp_t.state[c].j_2, J_2) != 1)
                    continue;
                  if (sp_t.state[a].l == 0 && sp_t.state[b].l == 0 && a == b &&
                      b == c && std::abs(T_2z) == 3)
                    continue;
                  state_t.a = a;
                  state_t.b = b;
                  state_t.c = c;
                  state_t.Jab_2 = J2ab;
                  state_t.Tab_2 = T2ab;
                  state_t.J_2 = J_2;
                  state_t.T_2 = T_2;
                  state_t.T_2z = T_2z;
                  state_t.E =
                      sp_t.state[a].E + sp_t.state[b].E + sp_t.state[c].E;
                  state_t.Par = this->phase(state_t.E);
                  state[chan_num].push_back(state_t);
                }
              }
            }
          }
        }
      }
    }
  }

  for (int i = 0; i < chan_size; i++) {
    stable_sort(state[i].begin(), state[i].end(), [](State_TP a, State_TP b) {
      return (1000000 * a.E + 100000 * a.Tab_2 + 10000 * a.Jab_2 + 100 * a.a +
              a.b + 0.01 * a.c + 0.001 * a.T_2z) <
             (1000000 * b.E + 100000 * b.Tab_2 + 10000 * b.Jab_2 + 100 * b.a +
              b.b + 0.01 * b.c + 0.001 * b.T_2z);
    });
    int index = 0;
    for (int j = 0; j < state[i].size(); j++) {
      state[i][j].index = index;
      index++;
    }
  }

  channel.resize(chan_size);

  for (int i = 0; i < chan_size; i++) {
    if (state[i].size() == 0) {
      channel[i].num = 0;
      channel[i].J_2 = 0;
      channel[i].T_2 = 0;
      channel[i].E_min = 0;
      channel[i].E_max = 0;
    } else {
      channel[i].num = state[i].size();
      channel[i].J_2 = state[i][0].J_2;
      channel[i].T_2 = state[i][0].T_2;
      int E_min = state[i][0].E;
      int E_max = state[i][0].E;
      for (int j = 0; j < state[i].size(); j++) {
        if (E_min > state[i][j].E)
          E_min = state[i][j].E;
        if (E_max < state[i][j].E)
          E_max = state[i][j].E;
      }
      channel[i].E_min = E_min;
      channel[i].E_max = E_max;
    }
  }

  // mark the indexE
  for (int i = 0; i < chan_size; i++) {
    int Emin = channel[i].E_min;
    int Emax = channel[i].E_max;
    channel[i].indexE.resize(Emax + 1);
    channel[i].indexE_beg.resize(Emax + 1);
    int flag = 1;
    for (int j = 0; j < state[i].size(); j++) {
      channel[i].indexE[state[i][j].E] = state[i][j].index;
    }
    for (int E = Emin; E <= Emax; E++) {
      int indexE_beg = 0;
      for (int j = 0; j < state[i].size(); j++) {
        if (state[i][j].E == E) {
          indexE_beg = j;
          break;
        }
      }
      channel[i].indexE_beg[E] = indexE_beg;
    }
  }

  // calulate JT2channel
  JT2channel.resize(2); // 0: T=0.5; 1:T=1.5
  for (int i = 0; i <= 1; i++) {
    int T_2 = i * 2 + 1;
    JT2channel[i].resize(J2_max + 1);
    for (int j = 0; j < JT2channel[i].size(); j++)
      JT2channel[i][j] = -1;
    for (int J2 = 1; J2 <= J2_max; J2 = J2 + 2) {
      for (int k = 0; k < channel.size(); k++) {
        if (channel[k].T_2 == T_2 && channel[k].J_2 == J2) {
          JT2channel[i][J2] = k;
          break;
        }
      }
    }
  }
}

void TripleP::build_T(SingleP &sp_t, int a_beg, int a_end, int b_beg, int b_end,
                      int c_beg, int c_end) {
  int Ea = 2 * sp_t.state[a_end].n + sp_t.state[a_end].l;
  int Eb = 2 * sp_t.state[b_end].n + sp_t.state[b_end].l;
  int Ec = 2 * sp_t.state[c_end].n + sp_t.state[c_end].l;
  Emax = Ea + Eb + Ec;
  int j2a_max = sp_t.l_max * 2 + 1;
  J2_max = 3 * j2a_max;
  int sp_size = sp_t.state.size();
  int chan_size = J2_max + 1;
  state.resize(chan_size);

  State_TP state_t;

  int chan_num = 0;
  for (int T_2 = 1; T_2 <= 3; T_2 = T_2 + 2) {
    for (int J_2 = 1; J_2 <= J2_max; J_2 = J_2 + 2) {
      // chan_num++;
      int T2ab_min = 2;
      int T2ab_max = 2;
      if (T_2 == 1) {
        T2ab_min = 0;
      }
      int J2ab_min, J2ab_max;
      J2ab_min = std::max(0, J_2 - j2a_max);
      J2ab_max = J_2 + j2a_max;
      for (int T2ab = T2ab_min; T2ab <= T2ab_max; T2ab = T2ab + 2) {
        // for (int T_2z = -T_2; T_2z <= T_2; T_2z = T_2z + 2) {
        // cout<<T_2z<<endl;
        for (int J2ab = J2ab_min; J2ab <= J2ab_max; J2ab = J2ab + 2) {
          for (int a = a_beg; a <= a_end; a++) {
            // for(int a=0; a<1; a++){
            // if(sp_t.state[a].E>Emax)continue;
            for (int b = b_beg; b <= b_end; b++) {
              // for(int b=0; b<1; b++){
              // if(sp_t.state[a].E+sp_t.state[b].E>Emax)continue;
              if (Tri(J2ab, sp_t.state[a].j_2, sp_t.state[b].j_2) != 1)
                continue;
              for (int c = c_beg; c <= c_end; c++) {
                // for(int c=0; c<1; c++){
                // if(sp_t.state[a].E+sp_t.state[b].E+sp_t.state[c].E>Emax)continue;

                if (Tri(J2ab, sp_t.state[c].j_2, J_2) != 1)
                  continue;
                /*if (sp_t.state[a].l == 0 && sp_t.state[b].l == 0 && a == b
                  && b == c && std::abs(T_2z) == 3) continue;*/
                state_t.I = -1;
                state_t.F = -1;
                state_t.Ecut_Flag = -1;
                state_t.a = a;
                state_t.b = b;
                state_t.c = c;
                state_t.Jab_2 = J2ab;
                state_t.Tab_2 = T2ab;
                state_t.J_2 = J_2;
                state_t.T_2 = T_2;
                state_t.T_2z = 0;
                state_t.E = sp_t.state[a].E + sp_t.state[b].E + sp_t.state[c].E;
                state_t.Par = this->phase(state_t.E);
                state[chan_num].push_back(state_t);
              }
            }
            //}
          }
        }
      }
      if (state[chan_num].size() != 0)
        chan_num++;
    }
  }

  for (int i = 0; i < chan_size; i++) {
    stable_sort(state[i].begin(), state[i].end(), [](State_TP a, State_TP b) {
      return (1000000 * a.E + 100000 * a.Tab_2 + 10000 * a.Jab_2 + 100 * a.a +
              a.b + 0.01 * a.c + 0.001 * a.T_2z) <
             (1000000 * b.E + 100000 * b.Tab_2 + 10000 * b.Jab_2 + 100 * b.a +
              b.b + 0.01 * b.c + 0.001 * b.T_2z);
    });
    int index = 0;
    for (int j = 0; j < state[i].size(); j++) {
      state[i][j].index = index;
      index++;
    }
  }
  // cout<<"good--?"<<endl;
  auto it = state.begin();
  while (it != state.end()) {
    if (it->size() == 0) {
      it = state.erase(it);
    } else {
      it++;
    }
  }
  chan_size = state.size();

  int all_size_T_t = 0;
  for (int i = 0; i < chan_size; i++) {
    all_size_T_t += state[i].size();
  }
  this->all_size_T = all_size_T_t;
  /*
    channel.resize(chan_size);

    for (int i = 0; i < chan_size; i++) {
      if (state[i].size() == 0) {
        channel[i].num = 0;
        channel[i].J_2 = 0;
        channel[i].T_2 = 0;
        channel[i].E_min = 0;
        channel[i].E_max = 0;
        // channel.erase(i);
      } else {
        channel[i].num = state[i].size();
        channel[i].J_2 = state[i][0].J_2;
        channel[i].T_2 = state[i][0].T_2;
        int E_min = state[i][0].E;
        int E_max = state[i][0].E;
        for (int j = 0; j < state[i].size(); j++) {
          if (E_min > state[i][j].E)
            E_min = state[i][j].E;
          if (E_max < state[i][j].E)
            E_max = state[i][j].E;
        }
        channel[i].E_min = E_min;
        channel[i].E_max = E_max;
      }
    }

    // mark the indexE
    for (int i = 0; i < chan_size; i++) {
      int Emin = channel[i].E_min;
      int Emax = channel[i].E_max;
      channel[i].indexE.resize(Emax + 1);
      channel[i].indexE_beg.resize(Emax + 1);
      int flag = 1;
      for (int j = 0; j < state[i].size(); j++) {
        channel[i].indexE[state[i][j].E] = state[i][j].index;
      }
      for (int E = Emin; E <= Emax; E++) {
        int indexE_beg = 0;
        for (int j = 0; j < state[i].size(); j++) {
          if (state[i][j].E == E) {
            indexE_beg = j;
            break;
          }
        }
        channel[i].indexE_beg[E] = indexE_beg;
      }
    }

    // calulate JT2channel
    JT2channel.resize(2); // 0: T=0.5; 1:T=1.5
    for (int i = 0; i <= 1; i++) {
      int T_2 = i * 2 + 1;
      JT2channel[i].resize(J2_max + 1);
      for (int j = 0; j < JT2channel[i].size(); j++)
        JT2channel[i][j] = -1;
      for (int J2 = 1; J2 <= J2_max; J2 = J2 + 2) {
        for (int k = 0; k < channel.size(); k++) {
          if (channel[k].T_2 == T_2 && channel[k].J_2 == J2) {
            JT2channel[i][J2] = k;
            break;
          }
        }
      }
    }

    // Find j2max_sp
    int j2max_sp_t = 1;
    int sp_min = std::min(a_beg, b_beg);
    sp_min = std::min(sp_min, c_beg);

    int sp_max = std::max(a_end, b_end);
    sp_max = std::max(sp_max, c_end);
    for (int i = sp_min; i <= sp_max; i++) {
      int j2 = sp_t.state[i].j_2;
      if (j2 > j2max_sp_t)
        j2max_sp_t = j2;
    }
    this->j2max_sp = j2max_sp_t;
    // Find J2_max
    int J2_max_t = 1;
    for (int i = 0; i < chan_size; i++) {
      if (channel[i].J_2 > J2_max_t)
        J2_max_t = channel[i].J_2;
    }
    this->J2_max_TP = J2_max_t;/**/
}

void TripleP::build_TTz(SingleP &sp_t, int a_beg, int a_end, int b_beg,
                        int b_end, int c_beg, int c_end) {
  int Ea = 2 * sp_t.state[a_end].n + sp_t.state[a_end].l;
  int Eb = 2 * sp_t.state[b_end].n + sp_t.state[b_end].l;
  int Ec = 2 * sp_t.state[c_end].n + sp_t.state[c_end].l;
  Emax = Ea + Eb + Ec;
  int j2a_max = sp_t.l_max * 2 + 1;
  J2_max = 3 * j2a_max;
  int sp_size = sp_t.state.size();
  int chan_size = J2_max + 1;
  state.resize(chan_size);

  State_TP state_t;
  int chan_num = 0;
  for (int T_2 = 1; T_2 <= 3; T_2 = T_2 + 2) {
    for (int J_2 = 1; J_2 <= J2_max; J_2 = J_2 + 2) {
      // chan_num++;
      int T2ab_min = 2;
      int T2ab_max = 2;
      if (T_2 == 1) {
        T2ab_min = 0;
      }
      int J2ab_min, J2ab_max;
      J2ab_min = std::max(0, J_2 - j2a_max);
      J2ab_max = J_2 + j2a_max;
      for (int T2ab = T2ab_min; T2ab <= T2ab_max; T2ab = T2ab + 2) {
        for (int T_2z = -T_2; T_2z <= T_2; T_2z = T_2z + 2) {
          // cout<<T_2z<<endl;
          for (int J2ab = J2ab_min; J2ab <= J2ab_max; J2ab = J2ab + 2) {
            for (int a = a_beg; a <= a_end; a++) {
              // for(int a=0; a<1; a++){
              // if(sp_t.state[a].E>Emax)continue;
              for (int b = b_beg; b <= b_end; b++) {
                // for(int b=0; b<1; b++){
                // if(sp_t.state[a].E+sp_t.state[b].E>Emax)continue;
                if (Tri(J2ab, sp_t.state[a].j_2, sp_t.state[b].j_2) != 1)
                  continue;
                for (int c = c_beg; c <= c_end; c++) {
                  // for(int c=0; c<1; c++){
                  // if(sp_t.state[a].E+sp_t.state[b].E+sp_t.state[c].E>Emax)continue;

                  if (Tri(J2ab, sp_t.state[c].j_2, J_2) != 1)
                    continue;
                  if (sp_t.state[a].l == 0 && sp_t.state[b].l == 0 && a == b &&
                      b == c && std::abs(T_2z) == 3)
                    continue;
                  state_t.a = a;
                  state_t.b = b;
                  state_t.c = c;
                  state_t.Jab_2 = J2ab;
                  state_t.Tab_2 = T2ab;
                  state_t.J_2 = J_2;
                  state_t.T_2 = T_2;
                  state_t.T_2z = T_2z;
                  state_t.E =
                      sp_t.state[a].E + sp_t.state[b].E + sp_t.state[c].E;
                  state_t.Par = this->phase(state_t.E);
                  state[chan_num].push_back(state_t);
                }
              }
            }
          }
        }
      }
      if (state[chan_num].size() != 0)
        chan_num++;
    }
  }

  for (int i = 0; i < chan_size; i++) {
    stable_sort(state[i].begin(), state[i].end(), [](State_TP a, State_TP b) {
      return (1000000 * a.E + 100000 * a.Tab_2 + 10000 * a.Jab_2 + 100 * a.a +
              a.b + 0.01 * a.c + 0.001 * a.T_2z) <
             (1000000 * b.E + 100000 * b.Tab_2 + 10000 * b.Jab_2 + 100 * b.a +
              b.b + 0.01 * b.c + 0.001 * b.T_2z);
    });
    int index = 0;
    for (int j = 0; j < state[i].size(); j++) {
      state[i][j].index = index;
      index++;
    }
  }
  // cout<<"good--?"<<endl;
  auto it = state.begin();
  while (it != state.end()) {
    if (it->size() == 0) {
      it = state.erase(it);
    } else {
      it++;
    }
  }
  chan_size = state.size();

  channel.resize(chan_size);

  for (int i = 0; i < chan_size; i++) {
    if (state[i].size() == 0) {
      channel[i].num = 0;
      channel[i].J_2 = 0;
      channel[i].T_2 = 0;
      channel[i].E_min = 0;
      channel[i].E_max = 0;
      // channel.erase(i);
    } else {
      channel[i].num = state[i].size();
      channel[i].J_2 = state[i][0].J_2;
      channel[i].T_2 = state[i][0].T_2;
      int E_min = state[i][0].E;
      int E_max = state[i][0].E;
      for (int j = 0; j < state[i].size(); j++) {
        if (E_min > state[i][j].E)
          E_min = state[i][j].E;
        if (E_max < state[i][j].E)
          E_max = state[i][j].E;
      }
      channel[i].E_min = E_min;
      channel[i].E_max = E_max;
    }
  }

  // mark the indexE
  for (int i = 0; i < chan_size; i++) {
    int Emin = channel[i].E_min;
    int Emax = channel[i].E_max;
    channel[i].indexE.resize(Emax + 1);
    channel[i].indexE_beg.resize(Emax + 1);
    int flag = 1;
    for (int j = 0; j < state[i].size(); j++) {
      channel[i].indexE[state[i][j].E] = state[i][j].index;
    }
    for (int E = Emin; E <= Emax; E++) {
      int indexE_beg = 0;
      for (int j = 0; j < state[i].size(); j++) {
        if (state[i][j].E == E) {
          indexE_beg = j;
          break;
        }
      }
      channel[i].indexE_beg[E] = indexE_beg;
    }
  }

  // calulate JT2channel
  JT2channel.resize(2); // 0: T=0.5; 1:T=1.5
  for (int i = 0; i <= 1; i++) {
    int T_2 = i * 2 + 1;
    JT2channel[i].resize(J2_max + 1);
    for (int j = 0; j < JT2channel[i].size(); j++)
      JT2channel[i][j] = -1;
    for (int J2 = 1; J2 <= J2_max; J2 = J2 + 2) {
      for (int k = 0; k < channel.size(); k++) {
        if (channel[k].T_2 == T_2 && channel[k].J_2 == J2) {
          JT2channel[i][J2] = k;
          break;
        }
      }
    }
  }

  // Find j2max_sp
  int j2max_sp_t = 1;
  int sp_min = std::min(a_beg, b_beg);
  sp_min = std::min(sp_min, c_beg);

  int sp_max = std::max(a_end, b_end);
  sp_max = std::max(sp_max, c_end);
  for (int i = sp_min; i <= sp_max; i++) {
    int j2 = sp_t.state[i].j_2;
    if (j2 > j2max_sp_t)
      j2max_sp_t = j2;
  }
  this->j2max_sp = j2max_sp_t;
}
void TripleP::Delet_Sort_T() {
  // call Select_FI_T first !!!
  // 1. Delete by .F = -1 && .I = -1
  // 2. Sort for index and channel
  //
  // Delete
  vector<vector<State_TP>> state_t;
  state_t = this->state;
  int chan_size = state_t.size();
  state.clear();
  state.resize(chan_size);

  for (int ch = 0; ch < chan_size; ch++) {
    int num_I = 0;
    int num_F = 0;
    for (int i = 0; i < state_t[ch].size(); i++) {
      if (state_t[ch][i].F > 0)
        num_F++;
      if (state_t[ch][i].I > 0)
        num_I++;
    }
    if (num_F > 0 && num_I > 0) {
      for (int i = 0; i < state_t[ch].size(); i++) {
        if (state_t[ch][i].F < 0 && state_t[ch][i].I < 0)
          continue;
        state[ch].push_back(state_t[ch][i]);
      }
    }
  }
  //
  // Sort
  // cout << "Good ? at TripleP::Delet_Sort_T" << endl;
  this->sort();
  // cout << "Good 2? at TripleP::Delet_Sort_T" << endl;
}
void TripleP::Delet_conf_T(SingleP &sp_t, vector<vector<int>> Sign,
                           vector<int> E_cut_gather,
                           vector<int> Same_Particles_gather) {
  // cout << "Start TripleP::Delet_conf_T " << endl;
  //
  // Delete
  int file_num = Sign.size();
  int size_ch = state.size();
  // this->E_cut_conf = E_cut;
  int TPconf_size = 0;

  // cout << "k = " << k << "\t E_cut_gather: " << E_cut_gather[k] << endl;
  for (int ch = 0; ch < size_ch; ch++) {
    int size = state[ch].size();
    for (int f = 0; f < size; f++) {
      if (state[ch][f].F < 0)
        continue;
      int ap = state[ch][f].a;
      int bp = state[ch][f].b;
      int cp = state[ch][f].c;
      int E_ap = sp_t.state[ap].E;
      int E_bp = sp_t.state[bp].E;
      int E_cp = sp_t.state[cp].E;
      for (int i = 0; i < state[ch].size(); i++) {
        if (state[ch][i].I < 0)
          continue;
        int a = state[ch][i].a;
        int b = state[ch][i].b;
        int c = state[ch][i].c;
        int E_a = sp_t.state[a].E;
        int E_b = sp_t.state[b].E;
        int E_c = sp_t.state[c].E;
        int Etot;
        int flag = -1;

        int flag_2 = -1;
        for (int k = 0; k < file_num; k++) {
          switch (Same_Particles_gather[k]) {
          case -1:
            flag_2 = 1;
            break;
          case 1:
            if (ap == a)
              flag_2 = 1;
            break;
          case 2:
            if (bp == b)
              flag_2 = 1;
            break;
          case 3:
            if (cp == c)
              flag_2 = 1;
            break;
          case 12:
            if (ap == a && bp == b)
              flag_2 = 1;
            break;
          case 13:
            if (ap == a && cp == c)
              flag_2 = 1;
            break;
          case 23:
            if (bp == b && cp == c)
              flag_2 = 1;
            break;
          case 123:
            if (ap == a && bp == b && cp == c)
              flag_2 = 1;
            break;
          default:
            cout << "Wrong happened here @ TripleP::build_conf_T "
                    "Same_Particles_gather[k] = "
                 << k << "\t" << Same_Particles_gather[k] << endl;
            exit(0);
          }
          if (flag_2 < 0)
            continue;
          if (this->FI_compare(state[ch][f].F, state[ch][i].I, k) < 0)
            continue;
          Etot = E_ap * Sign[k][0] + E_bp * Sign[k][1] + E_cp * Sign[k][2];
          Etot += E_a * Sign[k][3] + E_b * Sign[k][4] + E_c * Sign[k][5];
          if (std::abs(Etot) <= E_cut_gather[k]) {
            flag = 1;
          }
        }

        if (flag > 0) {
          state[ch][i].Ecut_Flag = 1;
          state[ch][f].Ecut_Flag = 1;
        }
      }
    }
  }

  vector<vector<State_TP>> state_t;
  state_t = this->state;
  int chan_size = state_t.size();
  state.clear();
  state.resize(chan_size);

  for (int ch = 0; ch < chan_size; ch++) {
    for (int i = 0; i < state_t[ch].size(); i++) {
      if (state_t[ch][i].Ecut_Flag != 1)
        continue;
      state[ch].push_back(state_t[ch][i]);
    }
  }

  auto it = state.begin();
  while (it != state.end()) {
    if (it->size() == 0) {
      it = state.erase(it);
    } else {
      it++;
    }
  }
  this->sort();
}
void TripleP::Delet_Sort_Tz() {
  // call Select_FI_T first !!!
  // 1. Delete by .F = -1 && .I = -1
  // 2. Sort for index and channel
  //
  // Delete
  vector<vector<State_TP>> state_Tz_t;
  state_Tz_t = this->state_Tz;
  int chan_size = state_Tz_t.size();
  state_Tz.clear();
  state_Tz.resize(chan_size);

  for (int ch = 0; ch < chan_size; ch++) {
    int num_I = 0;
    int num_F = 0;
    for (int i = 0; i < state_Tz_t[ch].size(); i++) {
      if (state_Tz_t[ch][i].F > 0)
        num_F++;
      if (state_Tz_t[ch][i].I > 0)
        num_I++;
    }
    if (num_F > 0 && num_I > 0) {
      for (int i = 0; i < state_Tz_t[ch].size(); i++) {
        if (state_Tz_t[ch][i].F < 0 && state_Tz_t[ch][i].I < 0)
          continue;
        state_Tz[ch].push_back(state_Tz_t[ch][i]);
      }
    }
  }
  this->sort_Tz();
}
void TripleP::Delet_conf_Tz(SingleP &sp_t, vector<vector<int>> Sign,
                            vector<int> E_cut_gather,
                            vector<int> Same_Particles_gather) {
  // cout << "Start TripleP::Delet_conf_Tz " << endl;

  // call build_conf_T first !!!
  // 1. Delete by .E < TPconf_E[ch].E_I && TPconf_E[ch].E_I
  // 2. Sort for index and channel
  //
  // Delete
  int file_num = Sign.size();
  int ch_num = state_Tz.size();
  /*cout << "Start TripleP::build_conf_Tz_S file_num : " << file_num
       << "\t ch_num :" << ch_num << endl;/**/

  int TPconf_size = 0;

  for (int ch = 0; ch < ch_num; ch++) {
    int size = state_Tz[ch].size();
    for (int f = 0; f < size; f++) {
      if (state_Tz[ch][f].F < 0)
        continue;
      int ap = state_Tz[ch][f].a;
      int bp = state_Tz[ch][f].b;
      int cp = state_Tz[ch][f].c;
      int E_ap = sp_t.state_Tz[ap].E;
      int E_bp = sp_t.state_Tz[bp].E;
      int E_cp = sp_t.state_Tz[cp].E;
      for (int i = 0; i < size; i++) {
        if (state_Tz[ch][i].I < 0)
          continue;
        int a = state_Tz[ch][i].a;
        int b = state_Tz[ch][i].b;
        int c = state_Tz[ch][i].c;
        int E_a = sp_t.state_Tz[a].E;
        int E_b = sp_t.state_Tz[b].E;
        int E_c = sp_t.state_Tz[c].E;

        int flag = -1;
        for (int k = 0; k < file_num; k++) {
          int flag_2 = -1;
          switch (Same_Particles_gather[k]) {
          case -1:
            flag_2 = 1;
            break;
          case 1:
            if (ap == a)
              flag_2 = 1;
            break;
          case 2:
            if (bp == b)
              flag_2 = 1;
            break;
          case 3:
            if (cp == c)
              flag_2 = 1;
            break;
          case 12:
            if (ap == a && bp == b)
              flag_2 = 1;
            break;
          case 13:
            if (ap == a && cp == c)
              flag_2 = 1;
            break;
          case 23:
            if (bp == b && cp == c)
              flag_2 = 1;
            break;
          case 123:
            if (ap == a && bp == b && cp == c)
              flag_2 = 1;
            break;
          default:
            cout << "Wrong happened here @ TripleP::build_conf_Tz "
                    "Same_Particles_gather[k] = "
                 << k << "\t" << Same_Particles_gather[k] << endl;
            exit(0);
          }
          if (flag_2 < 0)
            continue;
          if (this->FI_compare(state_Tz[ch][f].F, state_Tz[ch][i].I, k) < 0)
            continue;
          int Etot = E_ap * Sign[k][0] + E_bp * Sign[k][1] + E_cp * Sign[k][2];
          Etot += E_a * Sign[k][3] + E_b * Sign[k][4] + E_c * Sign[k][5];

          if (std::abs(Etot) <= E_cut_gather[k]) {
            flag = 1;
          }
        }
        if (flag > 0) {
          state_Tz[ch][i].Ecut_Flag = 1;
          state_Tz[ch][f].Ecut_Flag = 1;
        }
      }
    }
  }

  vector<vector<State_TP>> state_Tz_t;
  state_Tz_t = this->state_Tz;
  int chan_size = state_Tz.size();
  state_Tz.clear();
  state_Tz.resize(chan_size);
  // cout << "?? good chan_size = " << chan_size << endl;
  for (int ch = 0; ch < chan_size; ch++) {
    int E_FI_max = 0;

    for (int i = 0; i < state_Tz_t[ch].size(); i++) {
      if (state_Tz_t[ch][i].Ecut_Flag != 1)
        continue;
      state_Tz[ch].push_back(state_Tz_t[ch][i]);
    }
  }
  auto it = state_Tz.begin();
  while (it != state_Tz.end()) {
    if (it->size() == 0) {
      it = state_Tz.erase(it);
    } else {
      it++;
    }
  }
  // cout << "?? good " << endl;
  this->sort_Tz();
}
void TripleP::Select_FI_T(SingleP &sp_t,
                          vector<vector<vector<int>>> allow_FI_gather) {
  // orbits which particles can run :
  //  min: allow_FI[0][a', b', c', a, b, c]
  //  max: allow_FI[1][a', b', c', a, b, c]
  int file_num = allow_FI_gather.size();
  int sp_min_t = 100;
  int sp_max_t = -1;

  int Flag = 0;
  for (int k = 0; k < file_num; k++) {
    Flag = std::pow(10, k);

    vector<vector<int>> allow_FI;
    allow_FI = allow_FI_gather[k];
    // cout << "\ta'\t b'\t c'\t a\t b\t c\t" << endl;
    /*for (int i = 0; i < 2; i++) {
      if (i == 0)
        cout << "min: \t";
      if (i == 1)
        cout << "max: \t";
      for (int j = 0; j < 6; j++) {
        cout << allow_FI[i][j] << "\t";
        if (allow_FI[i][j] > sp_max_t)
          sp_max_t = allow_FI[i][j];
        if (allow_FI[i][j] < sp_min_t)
          sp_min_t = allow_FI[i][j];
      }
      cout << endl;
    }/**/

    int ap_min = allow_FI[0][0];
    int bp_min = allow_FI[0][1];
    int cp_min = allow_FI[0][2];
    int a_min = allow_FI[0][3];
    int b_min = allow_FI[0][4];
    int c_min = allow_FI[0][5];

    int ap_max = allow_FI[1][0];
    int bp_max = allow_FI[1][1];
    int cp_max = allow_FI[1][2];
    int a_max = allow_FI[1][3];
    int b_max = allow_FI[1][4];
    int c_max = allow_FI[1][5];

    for (int ch = 0; ch < state.size(); ch++) {
      for (int i = 0; i < state[ch].size(); i++) {
        int a = state[ch][i].a;
        int b = state[ch][i].b;
        int c = state[ch][i].c;
        if (a < ap_min || a > ap_max || b < bp_min || b > bp_max ||
            c < cp_min || c > cp_max) {
          int unuse;
        } else {
          if (state[ch][i].F < 0) {
            state[ch][i].F = Flag;
          } else {
            state[ch][i].F += Flag;
          }
        }
        if (a < a_min || a > a_max || b < b_min || b > b_max || c < c_min ||
            c > c_max) {
          int unuse;
        } else {
          if (state[ch][i].I < 0) {
            state[ch][i].I = Flag;
          } else {
            state[ch][i].I += Flag;
          }
        }
      }
    }
  }
  // Find j2max_sp
  int j2max_sp_t = 1;
  int sp_min = sp_min_t;
  int sp_max = sp_max_t;

  for (int i = sp_min; i <= sp_max; i++) {
    int j2 = sp_t.state[i].j_2;
    if (j2 > j2max_sp_t)
      j2max_sp_t = j2;
  }
  this->j2max_sp = j2max_sp_t;
}

void TripleP::Select_FI_Tz(vector<vector<vector<int>>> allow_FI_gather) {
  // orbits which particles can run :
  //  min: allow_FI[0][a', b', c', a, b, c]
  //  max: allow_FI[1][a, b, c, a, b, c]
  int file_num = allow_FI_gather.size();
  int Flag = 0;
  for (int k = 0; k < file_num; k++) {
    Flag = std::pow(10, k);
    vector<vector<int>> allow_FI;
    allow_FI = allow_FI_gather[k];
    int ap_min = allow_FI[0][0] * 2;
    int bp_min = allow_FI[0][1] * 2;
    int cp_min = allow_FI[0][2] * 2;
    int a_min = allow_FI[0][3] * 2;
    int b_min = allow_FI[0][4] * 2;
    int c_min = allow_FI[0][5] * 2;

    int ap_max = allow_FI[1][0] * 2 + 1;
    int bp_max = allow_FI[1][1] * 2 + 1;
    int cp_max = allow_FI[1][2] * 2 + 1;
    int a_max = allow_FI[1][3] * 2 + 1;
    int b_max = allow_FI[1][4] * 2 + 1;
    int c_max = allow_FI[1][5] * 2 + 1;

    for (int ch = 0; ch < state_Tz.size(); ch++) {
      for (int i = 0; i < state_Tz[ch].size(); i++) {
        int a = state_Tz[ch][i].a;
        int b = state_Tz[ch][i].b;
        int c = state_Tz[ch][i].c;
        if (a < ap_min || a > ap_max || b < bp_min || b > bp_max ||
            c < cp_min || c > cp_max) {
          int unuse;
        } else {
          if (state_Tz[ch][i].F < 0) {
            state_Tz[ch][i].F = Flag;
          } else {
            state_Tz[ch][i].F += Flag;
          }
        }
        if (a < a_min || a > a_max || b < b_min || b > b_max || c < c_min ||
            c > c_max) {
          int unuse;
        } else {
          if (state_Tz[ch][i].I < 0) {
            state_Tz[ch][i].I = Flag;
          } else {
            state_Tz[ch][i].I += Flag;
          }
        }
      }
    }
  }
}

int TripleP::FI_compare(int F, int I, int k) {
  // return 0: different; n>0 : allowed configure n times
  int n = -1;
  std::string F_s = std::to_string(F);
  std::string I_s = std::to_string(I);
  int size_F = F_s.size();
  int size_I = I_s.size();
  int size = std::min(size_F, size_I);
  if (k >= size)
    return n;
  /*if (k >= size) {
    cout << "Wrong happend at TripleP::FI_compare" << endl;
    cout << "F_s : " << F_s << "\t I_s : " << I_s << endl;
    cout << "F_s.length() = " << F_s.length() << "\t size = " << size
         << "\t k = " << k << endl;
    exit(0);
  }/**/
  int k_I = size_I - 1 - k;
  int k_F = size_F - 1 - k;
  char I_c = I_s[k_I];
  if ((F_s[k_F] == I_s[k_I]) && I_c != '0') {
    n = 1;
    // cout<<endl;
    // cout<<"k_F : "<<k_F<<"\t F_s : "<<F_s[k_F]<<"\t k_I:"<<k_I<<"\t
    // I_s[k_I]:"<<I_s[k_I]<<endl;
  }

  /*if ( size > 1 && n ==-1)
    cout << "F_s : " << F_s << "\t I_s : " << I_s << "k : " << k
         << "\t n= " << n << endl;/**/
  return n;
}

int TripleP::FI_compare(int F, int I) {
  // return 0: different; n>0 : allowed configure n times
  int n = -1;
  std::string F_s = std::to_string(F);
  std::string I_s = std::to_string(I);
  int size_F = F_s.size();
  int size_I = I_s.size();
  int size = std::min(size_F, size_I);

  /*if (k >= size) {
    cout << "Wrong happend at TripleP::FI_compare" << endl;
    cout << "F_s : " << F_s << "\t I_s : " << I_s << endl;
    cout << "F_s.length() = " << F_s.length() << "\t size = " << size
         << "\t k = " << k << endl;
    exit(0);
  }/**/
  for (int k = 0; k < size; k++) {
    int k_I = size_I - 1 - k;
    int k_F = size_F - 1 - k;
    char I_c = I_s[k_I];
    if ((F_s[k_F] == I_s[k_I]) && I_c != '0') {
      n = 1;
      break;
      // cout<<endl;
      // cout<<"k_F : "<<k_F<<"\t F_s : "<<F_s[k_F]<<"\t k_I:"<<k_I<<"\t
      // I_s[k_I]:"<<I_s[k_I]<<endl;
    }
  }

  /*if ( size > 1 && n ==-1)
    cout << "F_s : " << F_s << "\t I_s : " << I_s << "k : " << k
         << "\t n= " << n << endl;/**/
  return n;
}

void TripleP::build_conf_T(SingleP &sp_t, vector<vector<int>> Sign,
                           vector<int> E_cut_gather,
                           vector<int> Same_Particles_gather) {
  // build TPconf
  // | E_{omega} + E_{hole} - E_{particle} | < Ecut
  // E_Sign [a', b', c', a, b, c] = -1 {particle} +1 {omega or hole}
  /*int s_ap = E_Sign[0];
  int s_bp = E_Sign[1];
  int s_cp = E_Sign[2];
  int s_a = E_Sign[3];
  int s_b = E_Sign[4];
  int s_c = E_Sign[5];*/

  /*int A_t = 11;
  int B_t = 1001;
  cout<<"A_t : "<<A_t<<"\t B_t"<<B_t<<endl;
  for (int k = 0; k < 4; k++) {
    cout<<"\t\t \t k ="<<k<<"\t compare : "<<this->FI_compare(A_t, B_t,
  k)<<endl;
  }
  //exit(0);/**/
  int file_num = Sign.size();

  TPconf.resize(state.size());
  TPconf_k.resize(file_num);
  // this->E_cut_conf = E_cut;
  int TPconf_size = 0;
  for (int k = 0; k < file_num; k++) {
    TPconf_k[k].resize(state.size());
    // cout << "k = " << k << "\t E_cut_gather: " << E_cut_gather[k] << endl;
    for (int ch = 0; ch < state.size(); ch++) {
      for (int f = 0; f < state[ch].size(); f++) {
        if (state[ch][f].F < 0)
          continue;
        int ap = state[ch][f].a;
        int bp = state[ch][f].b;
        int cp = state[ch][f].c;
        int E_ap = sp_t.state[ap].E;
        int E_bp = sp_t.state[bp].E;
        int E_cp = sp_t.state[cp].E;
        for (int i = 0; i < state[ch].size(); i++) {
          if (state[ch][i].I < 0)
            continue;
          int a = state[ch][i].a;
          int b = state[ch][i].b;
          int c = state[ch][i].c;
          int E_a = sp_t.state[a].E;
          int E_b = sp_t.state[b].E;
          int E_c = sp_t.state[c].E;
          int Etot;
          int flag = -1;

          int flag_2 = -1;
          switch (Same_Particles_gather[k]) {
          case -1:
            flag_2 = 1;
            break;
          case 1:
            if (ap == a)
              flag_2 = 1;
            break;
          case 2:
            if (bp == b)
              flag_2 = 1;
            break;
          case 3:
            if (cp == c)
              flag_2 = 1;
            break;
          case 12:
            if (ap == a && bp == b)
              flag_2 = 1;
            break;
          case 13:
            if (ap == a && cp == c)
              flag_2 = 1;
            break;
          case 23:
            if (bp == b && cp == c)
              flag_2 = 1;
            break;
          case 123:
            if (ap == a && bp == b && cp == c)
              flag_2 = 1;
            break;
          default:
            cout << "Wrong happened here @ TripleP::build_conf_T "
                    "Same_Particles_gather[k] = "
                 << k << "\t" << Same_Particles_gather[k] << endl;
            exit(0);
          }
          if (flag_2 < 0)
            continue;

          if (this->FI_compare(state[ch][f].F, state[ch][i].I, k) < 0)
            continue;
          Etot = E_ap * Sign[k][0] + E_bp * Sign[k][1] + E_cp * Sign[k][2];
          Etot += E_a * Sign[k][3] + E_b * Sign[k][4] + E_c * Sign[k][5];
          if (std::abs(Etot) <= E_cut_gather[k]) {
            flag = 1;
          }

          if (flag > 0) {
            state[ch][i].Ecut_Flag = 1;
            state[ch][f].Ecut_Flag = 1;
            State_Conf conf_t;
            conf_t.F = f;
            conf_t.I = i;
            conf_t.E_ex = -1;
            TPconf[ch].push_back(conf_t);
            TPconf_k[k][ch].push_back(conf_t);
            TPconf_size++;
          }
        }
      }
    }
  }
  /*
  std::ofstream file;
  file.open("data/TPconf.dat");
  file << "Total size is : " << TPconf_size << endl;
  file << "File num : " << TPconf_k.size() << endl;
  for (int k = 0; k < TPconf_k.size(); k++) {
    file << "\t === file No. " << k << " ===" << endl;
    for (int ch = 0; ch < TPconf_k[k].size(); ch++) {
      file << "\t \t ch : " << ch << "\t size : " << TPconf_k[k][ch].size()
           << endl;
    }
  }
  for (int ch = 0; ch < TPconf.size(); ch++) {
    file << "\t == channel : " << ch << "\t == " << endl;
    for (int i = 0; i < TPconf[ch].size(); i++) {
      file << i << "\t" << TPconf[ch][i].F << "\t" << TPconf[ch][i].I << endl;
    }
  }
  file.close();/**/

  // &  build TPconf_E
  TPconf_E.resize(file_num);
  for (int k = 0; k < file_num; k++) {
    TPconf_E[k].resize(state.size());
    for (int ch = 0; ch < state.size(); ch++) {
      if (TPconf_k[k][ch].size() < 1)
        continue;

      for (int i = 0; i < TPconf_k[k][ch].size(); i++) {
        int F = TPconf_k[k][ch][i].F;
        int I = TPconf_k[k][ch][i].I;
        if (this->FI_compare(state[ch][F].F, state[ch][I].I, k) < 0)
          continue;
        int E_f = state[ch][F].E;
        int E_i = state[ch][I].E;
        State_TPconf_E TPconf_E_t;
        TPconf_E_t.E_F = E_f;
        TPconf_E_t.E_I = E_i;
        if (i == 0)
          TPconf_E[k][ch].push_back(TPconf_E_t);
        int flag = 1;
        for (int j = 0; j < TPconf_E[k][ch].size(); j++) {
          if (TPconf_E_t == TPconf_E[k][ch][j]) {
            flag = 0;
            break;
          }
        }
        if (flag == 1)
          TPconf_E[k][ch].push_back(TPconf_E_t);
      }
    }
  }

  /*file.open("data/TPconf_E.dat");
  file << "TPconf_E : " << endl;
  for (int k = 0; k < file_num; k++) {
    for (int ch = 0; ch < TPconf_E[k].size(); ch++) {
      file << "\t ---- File_num : " << k << " ----- ";
      file << "\t == channel : " << ch << "\t == " << endl;
      for (int i = 0; i < TPconf_E[k][ch].size(); i++) {
        file << i << "\t" << TPconf_E[k][ch][i].E_F << "\t"
             << TPconf_E[k][ch][i].E_I << endl;
      }
    }
  }
  file.close();/**/
  // &  build TPconf_E_tot
  TPconf_E_tot.resize(file_num);
  for (int k = 0; k < file_num; k++) {
    for (int ch = 0; ch < state.size(); ch++) {
      if (TPconf_E[k][ch].size() < 1)
        continue;
      if (TPconf_E_tot[k].size() == 0) {
        TPconf_E_tot[k].push_back(TPconf_E[k][ch][0]);
      }
      for (int i = 0; i < TPconf_E[k][ch].size(); i++) {
        int flag = 1;
        for (int j = 0; j < TPconf_E_tot[k].size(); j++) {
          if (TPconf_E[k][ch][i] == TPconf_E_tot[k][j]) {
            flag = 0;
            break;
          }
        }
        if (flag == 1)
          TPconf_E_tot[k].push_back(TPconf_E[k][ch][i]);
      }
    }
  }
  /*file.open("data/TPconf_E_tot.dat");
  file << "TPconf_E_tot : " << endl;
  for (int k = 0; k < file_num; k++) {
    file << "\t ---- File_num : " << k << " ----- " << endl;
    for (int i = 0; i < TPconf_E_tot[k].size(); i++) {
      file << i << "\t" << TPconf_E_tot[k][i].E_F << "\t"
           << TPconf_E_tot[k][i].E_I << endl;
    }
  }
  file.close();/**/
}

void TripleP::build_conf_T_S(SingleP &sp_t, vector<vector<int>> Sign,
                             vector<int> E_cut_gather,
                             vector<int> Same_Particles_gather) {
  // build TPconf
  // | E_{omega} + E_{hole} - E_{particle} | < Ecut
  // E_Sign [a', b', c', a, b, c] = -1 {particle} +1 {omega or hole}
  /*int s_ap = E_Sign[0];
  int s_bp = E_Sign[1];
  int s_cp = E_Sign[2];
  int s_a = E_Sign[3];
  int s_b = E_Sign[4];
  int s_c = E_Sign[5];*/

  int file_num = Sign.size();
  int size_ch = state.size();
  TPconf_S.cal_init(size_ch, 0);
  // this->E_cut_conf = E_cut;
  int TPconf_size = 0;

  // cout << "k = " << k << "\t E_cut_gather: " << E_cut_gather[k] << endl;
  for (int ch = 0; ch < size_ch; ch++) {
    int size = state[ch].size();
    TPconf_S.cal_beg(size, size);
    for (int f = 0; f < size; f++) {
      TPconf_S.cal_row_off(ch, f);
      if (state[ch][f].F < 0)
        continue;
      int ap = state[ch][f].a;
      int bp = state[ch][f].b;
      int cp = state[ch][f].c;
      int E_ap = sp_t.state[ap].E;
      int E_bp = sp_t.state[bp].E;
      int E_cp = sp_t.state[cp].E;
      for (int i = 0; i < state[ch].size(); i++) {
        if (state[ch][i].I < 0)
          continue;
        int a = state[ch][i].a;
        int b = state[ch][i].b;
        int c = state[ch][i].c;
        int E_a = sp_t.state[a].E;
        int E_b = sp_t.state[b].E;
        int E_c = sp_t.state[c].E;
        int Etot;
        int flag = -1;

        int flag_2 = -1;
        for (int k = 0; k < file_num; k++) {
          switch (Same_Particles_gather[k]) {
          case -1:
            flag_2 = 1;
            break;
          case 1:
            if (ap == a)
              flag_2 = 1;
            break;
          case 2:
            if (bp == b)
              flag_2 = 1;
            break;
          case 3:
            if (cp == c)
              flag_2 = 1;
            break;
          case 12:
            if (ap == a && bp == b)
              flag_2 = 1;
            break;
          case 13:
            if (ap == a && cp == c)
              flag_2 = 1;
            break;
          case 23:
            if (bp == b && cp == c)
              flag_2 = 1;
            break;
          case 123:
            if (ap == a && bp == b && cp == c)
              flag_2 = 1;
            break;
          default:
            cout << "Wrong happened here @ TripleP::build_conf_T "
                    "Same_Particles_gather[k] = "
                 << k << "\t" << Same_Particles_gather[k] << endl;
            exit(0);
          }
          if (flag_2 < 0)
            continue;
          if (this->FI_compare(state[ch][f].F, state[ch][i].I, k) < 0)
            continue;
          Etot = E_ap * Sign[k][0] + E_bp * Sign[k][1] + E_cp * Sign[k][2];
          Etot += E_a * Sign[k][3] + E_b * Sign[k][4] + E_c * Sign[k][5];
          if (std::abs(Etot) <= E_cut_gather[k]) {
            flag = 1;
          }
        }

        if (flag > 0) {
          state[ch][i].Ecut_Flag = 1;
          state[ch][f].Ecut_Flag = 1;
          TPconf_S.cal_col_ind(ch, f, i, 1);
        }
      }
    }
    TPconf_S.cal_end(ch);
  }

  /*
  std::ofstream file;
  file.open("data/TPconf.dat");
  file << "Total size is : " << TPconf_size << endl;
  file << "File num : " << TPconf_k.size() << endl;
  for (int k = 0; k < TPconf_k.size(); k++) {
    file << "\t === file No. " << k << " ===" << endl;
    for (int ch = 0; ch < TPconf_k[k].size(); ch++) {
      file << "\t \t ch : " << ch << "\t size : " << TPconf_k[k][ch].size()
           << endl;
    }
  }
  for (int ch = 0; ch < TPconf.size(); ch++) {
    file << "\t == channel : " << ch << "\t == " << endl;
    for (int i = 0; i < TPconf[ch].size(); i++) {
      file << i << "\t" << TPconf[ch][i].F << "\t" << TPconf[ch][i].I << endl;
    }
  }
  file.close();/**/

  // &  build TPconf_E
  TPconf_E.resize(file_num);
  for (int k = 0; k < file_num; k++) {
    TPconf_E[k].resize(state.size());
    for (int ch = 0; ch < state.size(); ch++) {
      if (TPconf_k[k][ch].size() < 1)
        continue;

      for (int i = 0; i < TPconf_k[k][ch].size(); i++) {
        int F = TPconf_k[k][ch][i].F;
        int I = TPconf_k[k][ch][i].I;
        if (this->FI_compare(state[ch][F].F, state[ch][I].I, k) < 0)
          continue;
        int E_f = state[ch][F].E;
        int E_i = state[ch][I].E;
        State_TPconf_E TPconf_E_t;
        TPconf_E_t.E_F = E_f;
        TPconf_E_t.E_I = E_i;
        if (i == 0)
          TPconf_E[k][ch].push_back(TPconf_E_t);
        int flag = 1;
        for (int j = 0; j < TPconf_E[k][ch].size(); j++) {
          if (TPconf_E_t == TPconf_E[k][ch][j]) {
            flag = 0;
            break;
          }
        }
        if (flag == 1)
          TPconf_E[k][ch].push_back(TPconf_E_t);
      }
    }
  }

  /*file.open("data/TPconf_E.dat");
  file << "TPconf_E : " << endl;
  for (int k = 0; k < file_num; k++) {
    for (int ch = 0; ch < TPconf_E[k].size(); ch++) {
      file << "\t ---- File_num : " << k << " ----- ";
      file << "\t == channel : " << ch << "\t == " << endl;
      for (int i = 0; i < TPconf_E[k][ch].size(); i++) {
        file << i << "\t" << TPconf_E[k][ch][i].E_F << "\t"
             << TPconf_E[k][ch][i].E_I << endl;
      }
    }
  }
  file.close();/**/
  // &  build TPconf_E_tot
  TPconf_E_tot.resize(file_num);
  for (int k = 0; k < file_num; k++) {
    for (int ch = 0; ch < state.size(); ch++) {
      if (TPconf_E[k][ch].size() < 1)
        continue;
      if (TPconf_E_tot[k].size() == 0) {
        TPconf_E_tot[k].push_back(TPconf_E[k][ch][0]);
      }
      for (int i = 0; i < TPconf_E[k][ch].size(); i++) {
        int flag = 1;
        for (int j = 0; j < TPconf_E_tot[k].size(); j++) {
          if (TPconf_E[k][ch][i] == TPconf_E_tot[k][j]) {
            flag = 0;
            break;
          }
        }
        if (flag == 1)
          TPconf_E_tot[k].push_back(TPconf_E[k][ch][i]);
      }
    }
  }
  /*file.open("data/TPconf_E_tot.dat");
  file << "TPconf_E_tot : " << endl;
  for (int k = 0; k < file_num; k++) {
    file << "\t ---- File_num : " << k << " ----- " << endl;
    for (int i = 0; i < TPconf_E_tot[k].size(); i++) {
      file << i << "\t" << TPconf_E_tot[k][i].E_F << "\t"
           << TPconf_E_tot[k][i].E_I << endl;
    }
  }
  file.close();/**/
}

void TripleP::build_conf_Tz(SingleP &sp_t, vector<vector<int>> Sign,
                            vector<int> E_cut_gather,
                            vector<int> Same_Particles_gather) {
  // build TPconf
  // | E_{omega} + E_{hole} - E_{particle} | < Ecut
  // E_Sign [a', b', c', a, b, c] = -1 {particle} +1 {omega or hole}
  /*int s_ap = E_Sign[0];
  int s_bp = E_Sign[1];
  int s_cp = E_Sign[2];
  int s_a = E_Sign[3];
  int s_b = E_Sign[4];
  int s_c = E_Sign[5];*/
  int file_num = Sign.size();
  int ch_num = state_Tz.size();
  /*cout << "Start TripleP::build_conf_Tz file_num : " << file_num
       << "\t ch_num :" << ch_num << endl;/**/
  TPconf_Tz.resize(ch_num);
  TPconf_k_Tz.resize(file_num);
  // this->E_cut_conf = E_cut;
  int TPconf_size = 0;
  for (int k = 0; k < file_num; k++) {
    TPconf_k_Tz[k].resize(ch_num);
    for (int ch = 0; ch < ch_num; ch++) {
      for (int f = 0; f < state_Tz[ch].size(); f++) {
        if (state_Tz[ch][f].F < 0)
          continue;
        int ap = state_Tz[ch][f].a;
        int bp = state_Tz[ch][f].b;
        int cp = state_Tz[ch][f].c;
        int E_ap = sp_t.state_Tz[ap].E;
        int E_bp = sp_t.state_Tz[bp].E;
        int E_cp = sp_t.state_Tz[cp].E;
        for (int i = 0; i < state_Tz[ch].size(); i++) {
          if (state_Tz[ch][i].I < 0)
            continue;
          int a = state_Tz[ch][i].a;
          int b = state_Tz[ch][i].b;
          int c = state_Tz[ch][i].c;
          int E_a = sp_t.state_Tz[a].E;
          int E_b = sp_t.state_Tz[b].E;
          int E_c = sp_t.state_Tz[c].E;
          // vector<int> Etot(file_num);

          int flag_2 = -1;
          switch (Same_Particles_gather[k]) {
          case -1:
            flag_2 = 1;
            break;
          case 1:
            if (ap == a)
              flag_2 = 1;
            break;
          case 2:
            if (bp == b)
              flag_2 = 1;
            break;
          case 3:
            if (cp == c)
              flag_2 = 1;
            break;
          case 12:
            if (ap == a && bp == b)
              flag_2 = 1;
            break;
          case 13:
            if (ap == a && cp == c)
              flag_2 = 1;
            break;
          case 23:
            if (bp == b && cp == c)
              flag_2 = 1;
            break;
          case 123:
            if (ap == a && bp == b && cp == c)
              flag_2 = 1;
            break;
          default:
            cout << "Wrong happened here @ TripleP::build_conf_Tz "
                    "Same_Particles_gather[k] = "
                 << k << "\t" << Same_Particles_gather[k] << endl;
            exit(0);
          }
          if (flag_2 < 0)
            continue;
          if (this->FI_compare(state_Tz[ch][f].F, state_Tz[ch][i].I, k) < 0)
            continue;
          int Etot = E_ap * Sign[k][0] + E_bp * Sign[k][1] + E_cp * Sign[k][2];
          Etot += E_a * Sign[k][3] + E_b * Sign[k][4] + E_c * Sign[k][5];
          int flag = -1;
          if (std::abs(Etot) <= E_cut_gather[k]) {
            state_Tz[ch][i].Ecut_Flag = 1;
            state_Tz[ch][f].Ecut_Flag = 1;
            State_Conf conf_t;
            conf_t.F = f;
            conf_t.I = i;
            conf_t.E_ex = -1;
            TPconf_Tz[ch].push_back(conf_t);
            TPconf_k_Tz[k][ch].push_back(conf_t);
            TPconf_size++;
          }
        }
      }
    }
  }
  /*
  for (int ch = 0; ch < ch_num; ch++) {
    cout << "ch : " << ch
         << "\t TPconf_Tz[ch].size() : " << TPconf_Tz[ch].size() << endl;
  }
  cout << "Finshed build_conf_Tz " << endl;/**/
  /*std::ofstream file;
  file.open("data/TPconf_Tz.dat");
  file << "Total_Tz size is : " << TPconf_size << endl;
  for (int ch = 0; ch < TPconf_Tz.size(); ch++) {
    file << "\t == channel : " << ch << "\t == " << endl;
    for (int i = 0; i < TPconf_Tz[ch].size(); i++) {
      file << i << "\t" << TPconf_Tz[ch][i].F << "\t" << TPconf_Tz[ch][i].I
           << endl;
    }
  }
  file.close();/**/
  // &  build TPconf_E
  TPconf_E_Tz.resize(file_num);
  for (int k = 0; k < file_num; k++) {
    TPconf_E_Tz[k].resize(state_Tz.size());
    for (int ch = 0; ch < state_Tz.size(); ch++) {
      if (TPconf_k_Tz[k][ch].size() < 1)
        continue;

      for (int i = 0; i < TPconf_k_Tz[k][ch].size(); i++) {
        int F = TPconf_k_Tz[k][ch][i].F;
        int I = TPconf_k_Tz[k][ch][i].I;
        if (this->FI_compare(state_Tz[ch][F].F, state_Tz[ch][I].I, k) < 0)
          continue;
        int E_f = state_Tz[ch][F].E;
        int E_i = state_Tz[ch][I].E;
        State_TPconf_E TPconf_E_t;
        TPconf_E_t.E_F = E_f;
        TPconf_E_t.E_I = E_i;
        if (i == 0)
          TPconf_E_Tz[k][ch].push_back(TPconf_E_t);
        int flag = 1;
        for (int j = 0; j < TPconf_E_Tz[k][ch].size(); j++) {
          if (TPconf_E_t == TPconf_E_Tz[k][ch][j]) {
            flag = 0;
            break;
          }
        }
        if (flag == 1)
          TPconf_E_Tz[k][ch].push_back(TPconf_E_t);
      }
    }
  }
  /*
    file.open("data/TPconf_E_Tz.dat");
    file << "TPconf_E_Tz : " << endl;
    for (int k = 0; k < file_num; k++) {
      for (int ch = 0; ch < TPconf_E_Tz[k].size(); ch++) {
        file << "\t ---- File_num : " << k << " ----- ";
        file << "\t == channel : " << ch << "\t == " << endl;
        for (int i = 0; i < TPconf_E_Tz[k][ch].size(); i++) {
          file << i << "\t" << TPconf_E_Tz[k][ch][i].E_F << "\t"
               << TPconf_E_Tz[k][ch][i].E_I << endl;
        }
      }
    }
    file.close();/**/
}

void TripleP::build_conf_Tz_S(SingleP &sp_t, vector<vector<int>> Sign,
                              vector<int> E_cut_gather,
                              vector<int> Same_Particles_gather) {
  // build TPconf
  // | E_{omega} + E_{hole} - E_{particle} | < Ecut
  // E_Sign [a', b', c', a, b, c] = -1 {particle} +1 {omega or hole}
  /*int s_ap = E_Sign[0];
  int s_bp = E_Sign[1];
  int s_cp = E_Sign[2];
  int s_a = E_Sign[3];
  int s_b = E_Sign[4];
  int s_c = E_Sign[5];*/
  int file_num = Sign.size();
  int ch_num = state_Tz.size();
  cout << "Start TripleP::build_conf_Tz_S file_num : " << file_num
       << "\t ch_num :" << ch_num << endl;
  TPconf_Tz_S.cal_init(ch_num, 0);

  int TPconf_size = 0;

  for (int ch = 0; ch < ch_num; ch++) {
    int size = state_Tz[ch].size();
    TPconf_Tz_S.cal_beg(size, size);
    for (int f = 0; f < size; f++) {
      TPconf_Tz_S.cal_row_off(ch, f);
      if (state_Tz[ch][f].F < 0)
        continue;
      int ap = state_Tz[ch][f].a;
      int bp = state_Tz[ch][f].b;
      int cp = state_Tz[ch][f].c;
      int E_ap = sp_t.state_Tz[ap].E;
      int E_bp = sp_t.state_Tz[bp].E;
      int E_cp = sp_t.state_Tz[cp].E;
      for (int i = 0; i < size; i++) {
        if (state_Tz[ch][i].I < 0)
          continue;
        int a = state_Tz[ch][i].a;
        int b = state_Tz[ch][i].b;
        int c = state_Tz[ch][i].c;
        int E_a = sp_t.state_Tz[a].E;
        int E_b = sp_t.state_Tz[b].E;
        int E_c = sp_t.state_Tz[c].E;

        int flag = -1;
        for (int k = 0; k < file_num; k++) {
          int flag_2 = -1;
          switch (Same_Particles_gather[k]) {
          case -1:
            flag_2 = 1;
            break;
          case 1:
            if (ap == a)
              flag_2 = 1;
            break;
          case 2:
            if (bp == b)
              flag_2 = 1;
            break;
          case 3:
            if (cp == c)
              flag_2 = 1;
            break;
          case 12:
            if (ap == a && bp == b)
              flag_2 = 1;
            break;
          case 13:
            if (ap == a && cp == c)
              flag_2 = 1;
            break;
          case 23:
            if (bp == b && cp == c)
              flag_2 = 1;
            break;
          case 123:
            if (ap == a && bp == b && cp == c)
              flag_2 = 1;
            break;
          default:
            cout << "Wrong happened here @ TripleP::build_conf_Tz "
                    "Same_Particles_gather[k] = "
                 << k << "\t" << Same_Particles_gather[k] << endl;
            exit(0);
          }
          if (flag_2 < 0)
            continue;
          if (this->FI_compare(state_Tz[ch][f].F, state_Tz[ch][i].I, k) < 0)
            continue;
          int Etot = E_ap * Sign[k][0] + E_bp * Sign[k][1] + E_cp * Sign[k][2];
          Etot += E_a * Sign[k][3] + E_b * Sign[k][4] + E_c * Sign[k][5];

          if (std::abs(Etot) <= E_cut_gather[k]) {
            flag = 1;
          }
        }
        if (flag > 0) {
          state_Tz[ch][i].Ecut_Flag = 1;
          state_Tz[ch][f].Ecut_Flag = 1;
          TPconf_Tz_S.cal_col_ind(ch, f, i, 1);
        }
      }
    }
    TPconf_Tz_S.cal_end(ch);
  }

  cout << "Finshed build_conf_Tz_S " << endl;
}

void TripleP::find_E_Jac_conf(vector<vector<int>> &E_cut_conf_t) {
  int file_num = this->TPconf_E_tot.size();
  if (file_num < 1) {
    cout << "Wrong @ TripleP::find_E_Jac_conf " << endl;
    cout << "\t this->TPconf_E_tot.size() <1 " << endl;
  }
  for (int k = 0; k < file_num; k++) {
    int size = TPconf_E_tot[k].size();
    {
      for (int i = 0; i < size; i++) {
        int E_F = TPconf_E_tot[k][i].E_F;
        int E_I = TPconf_E_tot[k][i].E_I;
        int flag = -1;
        int E_min = std::min(E_F, E_I);
        int E_Diff = 0;
        if (E_F <= E_I) {
          E_Diff = E_I - E_F;
          flag = 1;
        } else {
          E_Diff = E_F - E_I;
          flag = 2;
        }
        for (int E_1 = 0; E_1 <= E_min; E_1++) {
          int E_2 = E_1 + E_Diff;
          int E_F_t = -1;
          int E_I_t = -1;
          if (flag == 1) {
            E_F_t = E_1;
            E_I_t = E_2;
          } else {
            E_F_t = E_2;
            E_I_t = E_1;
          }
          int flag_2 = 1;
          int E_size = E_cut_conf_t.size();
          for (int j = 0; j < E_size; j++) {
            if ((E_cut_conf_t[j][0] == E_F_t) &&
                (E_cut_conf_t[j][1] == E_I_t)) {
              flag_2 = -1;
              break;
            }
          }
          // cout << "E_F_t : " << E_F_t << "\t E_I_t" << E_I_t
          //      << "\t flag_2 : " << flag_2 << endl;
          if (flag_2 == 1) {
            vector<int> vec_t(2);
            vec_t[0] = E_F_t;
            vec_t[1] = E_I_t;
            E_cut_conf_t.push_back(vec_t);
          }
        }
      }
    }
  }
  this->E_Jac_conf = E_cut_conf_t;
}

void TripleP::sort() {
  // Sort state
  auto it = state.begin();
  while (it != state.end()) {
    if (it->size() == 0) {
      it = state.erase(it);
    } else {
      it++;
    }
  }
  int chan_size = state.size();
  for (int i = 0; i < chan_size; i++) {
    stable_sort(state[i].begin(), state[i].end(), [](State_TP a, State_TP b) {
      return (1000000 * a.E + 100000 * a.Tab_2 + 10000 * a.Jab_2 + 100 * a.a +
              a.b + 0.01 * a.c + 0.001 * a.T_2z) <
             (1000000 * b.E + 100000 * b.Tab_2 + 10000 * b.Jab_2 + 100 * b.a +
              b.b + 0.01 * b.c + 0.001 * b.T_2z);
    });
    int index = 0;
    for (int j = 0; j < state[i].size(); j++) {
      state[i][j].index = index;
      index++;
    }
  }

  // cout<<"good? 1 "<<endl;
  // Rebuild Channel
  chan_size = state.size();
  channel.clear();

  int all_size_T_t = 0;
  for (int i = 0; i < chan_size; i++) {
    all_size_T_t += state[i].size();
  }
  this->all_size_T = all_size_T_t;

  channel.resize(chan_size);

  for (int i = 0; i < chan_size; i++) {
    if (state[i].size() == 0) {
      channel[i].num = 0;
      channel[i].J_2 = 0;
      channel[i].T_2 = 0;
      channel[i].E_min = 0;
      channel[i].E_max = 0;
      // channel.erase(i);
    } else {
      channel[i].num = state[i].size();
      channel[i].J_2 = state[i][0].J_2;
      channel[i].T_2 = state[i][0].T_2;
      int E_min = state[i][0].E;
      int E_max = state[i][0].E;
      for (int j = 0; j < state[i].size(); j++) {
        if (E_min > state[i][j].E)
          E_min = state[i][j].E;
        if (E_max < state[i][j].E)
          E_max = state[i][j].E;
      }
      channel[i].E_min = E_min;
      channel[i].E_max = E_max;
    }
  }
  // cout<<"good? 2 "<<endl;

  // mark the indexE
  for (int i = 0; i < chan_size; i++) {
    int Emin = channel[i].E_min;
    int Emax = channel[i].E_max;
    channel[i].indexE.resize(Emax + 1);
    channel[i].indexE_beg.resize(Emax + 1);
    int flag = 1;
    for (int j = 0; j < state[i].size(); j++) {
      channel[i].indexE[state[i][j].E] = state[i][j].index;
    }
    for (int E = Emin; E <= Emax; E++) {
      int indexE_beg = 0;
      for (int j = 0; j < state[i].size(); j++) {
        if (state[i][j].E == E) {
          indexE_beg = j;
          break;
        }
      }
      channel[i].indexE_beg[E] = indexE_beg;
    }
  }
  // mark channel_TP.E_min/max_F/I
  for (int ch = 0; ch < chan_size; ch++) {
    int Emin_f_t = 100;
    int Emin_i_t = 100;
    int Emax_f_t = 0;
    int Emax_i_t = 0;
    int Emax_t = 0;
    int Emin_t = 100;

    for (int i = 0; i < state[ch].size(); i++) {
      int E_t = state[ch][i].E;
      if (Emin_t > E_t)
        Emin_t = E_t;
      if (Emax_t < E_t)
        Emax_t = E_t;
      if (state[ch][i].F > 0) {
        Emin_f_t = Emin_t;
        Emax_f_t = Emax_t;
      }
      if (state[ch][i].I > 0) {
        Emin_i_t = Emin_t;
        Emax_i_t = Emax_t;
      }
    }

    this->channel[ch].E_min_F = Emin_f_t;
    this->channel[ch].E_min_I = Emin_i_t;
    this->channel[ch].E_max_F = Emax_f_t;
    this->channel[ch].E_max_I = Emax_i_t;
  }
  // cout<<"good? 3 "<<endl;

  // Find J2_max
  int J2_max_t = 1;
  for (int i = 0; i < chan_size; i++) {
    if (channel[i].J_2 > J2_max_t)
      J2_max_t = channel[i].J_2;
  }
  // cout<<"good? 4 "<<endl;

  this->J2_max_TP = J2_max_t;
  // cout<<"J2_max = "<<J2_max<<"\t J2_max_TP = "<<J2_max_TP<<endl;
  // calulate JT2channel
  JT2channel.resize(2); // 0: T=0.5; 1:T=1.5
  for (int i = 0; i <= 1; i++) {
    int T_2 = i * 2 + 1;
    JT2channel[i].resize(J2_max_TP + 1);
    for (int j = 0; j < JT2channel[i].size(); j++)
      JT2channel[i][j] = -1;
    for (int J2 = 1; J2 <= J2_max_TP; J2 = J2 + 2) {
      JT2channel[i][J2] = -1;
      for (int k = 0; k < channel.size(); k++) {
        if (channel[k].T_2 == T_2 && channel[k].J_2 == J2) {
          JT2channel[i][J2] = k;
          break;
        }
      }
    }
  }
  // cout<<"good? 5 "<<endl;

  int Emax_f = -1;
  int Emax_i = -1;
  int Emin_f = 100;
  int Emin_i = 100;
  for (int ch = 0; ch < state.size(); ch++) {
    for (int i = 0; i < state[ch].size(); i++) {
      int E = state[ch][i].E;
      if (state[ch][i].F > 0 && E > Emax_f) {
        Emax_f = E;
      }
      if (state[ch][i].F > 0 && E < Emin_f) {
        Emin_f = E;
      }
      if (state[ch][i].I > 0 && E > Emax_i) {
        Emax_i = E;
      }
      if (state[ch][i].I > 0 && E < Emin_i) {
        Emin_i = E;
      }
    }
  }
  if (Emax_f < 0 || Emax_i < 0 || Emin_f > 50 || Emin_i > 50) {
    cout << "Wrong happened at  TripleP::Select_FI_T " << endl;
    cout << "Emax : " << Emax_f << "\t " << Emax_i << "\t Emin : " << Emin_f
         << "\t" << Emin_i << endl;
    exit(0);
  }
  this->Emax_Final = Emax_f;
  this->Emax_Init = Emax_i;
  this->Emin_Final = Emin_f;
  this->Emin_Init = Emin_i;
  this->Emax = std::max(Emax_f, Emax_i);
  // cout<<"good? 6 "<<endl;
}

void TripleP::sort_Tz() {
  auto it = state_Tz.begin();
  while (it != state_Tz.end()) {
    if (it->size() == 0) {
      it = state_Tz.erase(it);
    } else {
      it++;
    }
  }
  int chan_size = state_Tz.size();
  for (int i = 0; i < chan_size; i++) {
    stable_sort(state_Tz[i].begin(), state_Tz[i].end(),
                [](State_TP a, State_TP b) {
                  return (100000 * a.E + 10000 * a.Jab_2 + 100 * a.a + a.b +
                          0.01 * a.c + 0.001 * a.Tab_2z + 0.0001 * a.T_2z) <
                         (100000 * b.E + 10000 * b.Jab_2 + 100 * b.a + b.b +
                          0.01 * b.c + 0.001 * b.Tab_2z + 0.0001 * b.T_2z);
                });
    int index = 0;
    for (int j = 0; j < state_Tz[i].size(); j++) {
      state_Tz[i][j].index = index;
      index++;
    }
  }

  // Rebuild channel
  chan_size = state_Tz.size();
  channel_Tz.clear();

  int all_size_Tz_t = 0;
  for (int i = 0; i < chan_size; i++) {
    all_size_Tz += chan_size;
  }
  this->all_size_Tz = all_size_Tz_t;

  channel_Tz.resize(chan_size);

  for (int i = 0; i < chan_size; i++) {
    if (state_Tz[i].size() == 0) {
      channel_Tz[i].num = 0;
      channel_Tz[i].J_2 = 0;
      channel_Tz[i].T_2z = 0;
      channel_Tz[i].E_min = 0;
      channel_Tz[i].E_max = 0;
    } else {
      channel_Tz[i].num = state_Tz[i].size();
      channel_Tz[i].J_2 = state_Tz[i][0].J_2;
      channel_Tz[i].T_2z = state_Tz[i][0].T_2z;
      int E_min = state_Tz[i][0].E;
      int E_max = state_Tz[i][0].E;
      for (int j = 0; j < state_Tz[i].size(); j++) {
        if (E_min > state_Tz[i][j].E)
          E_min = state_Tz[i][j].E;
        if (E_max < state_Tz[i][j].E)
          E_max = state_Tz[i][j].E;
      }
      channel_Tz[i].E_min = E_min;
      channel_Tz[i].E_max = E_max;
    }
  }

  // mark the indexE
  for (int i = 0; i < chan_size; i++) {
    int Emin = channel_Tz[i].E_min;
    int Emax = channel_Tz[i].E_max;
    channel_Tz[i].indexE.resize(Emax + 1);
    channel_Tz[i].indexE_beg.resize(Emax + 1);
    int flag = 1;
    for (int j = 0; j < state_Tz[i].size(); j++) {
      channel_Tz[i].indexE[state_Tz[i][j].E] = state_Tz[i][j].index;
    }
    for (int E = Emin; E <= Emax; E++) {
      int indexE_beg = 0;
      for (int j = 0; j < state_Tz[i].size(); j++) {
        if (state_Tz[i][j].E == E) {
          indexE_beg = j;
          break;
        }
      }
      channel_Tz[i].indexE_beg[E] = indexE_beg;
    }
  } /**/
}
int TripleP::FindC(int J2, int T2) {
  // Find Channel with J2 T2
  int Tx = 0;
  if (T2 == 3)
    Tx = 1; // 0: Tx=0.5; 1:Tx=1.5

  if (J2 > JT2channel[Tx].size()) {
    cout << "wrong at find channel in TripleP::FindC " << endl;
    cout << "\t J2 = " << J2 << "\t T2 = " << T2 << endl;
    exit(0);
  }
  return JT2channel[Tx][J2];
  /*if (size == 0)
    return -1;
  int front = 0;
  int end = size - 1;
  int mid = size / 2 - 1;
  if (T2 == 1) {
    end = mid;
    mid = floor((end - front) / 2) + front;
  } else {
    front = mid + 1;
    mid = floor((end - front) / 2) + front;
  }
  while (front < end && channel[mid].J_2 != J2) {
    if (channel[mid].J_2 < J2)
      front = mid + 1;
    if (channel[mid].J_2 > J2)
      end = mid - 1;
    mid = (front + end) / 2;
  }
  if (channel[mid].J_2 != J2) {
    cout << "wrong at find channel in TripleP::FindC " << endl;
    cout << "\t J2 = " << J2 << "\t T2 = " << T2 << endl;
    exit(0);
  } else {
    return mid;
  }/**/
}
int TripleP::FindC_Tz(int J2, int T_2z) {
  // Find Channel with J2 T2
  int res = -1;
  int size = channel_Tz.size();
  if (size == 0)
    return -1;
  for (int i = 0; i < size; i++) {
    if (J2 == channel_Tz[i].J_2 && T_2z == channel_Tz[i].T_2z) {
      res = i;
      break;
    }
  }
  if (res == -1) {
    cout << "wrong at find channel in TripleP::FindC_Tz " << endl;
    cout << "\t J2 = " << J2 << "\t T_2z = " << T_2z << endl;
    exit(0);
  }
  return res;
}

int TripleP::FindTP_Tz(int J2, int T_2z, int J2ab, int a, int b, int c) {
  int index = -1;
  int ch_Tz = FindC_Tz(J2, T_2z);
  int size = state_Tz[ch_Tz].size();
  for (int i = 0; i < size; i++) {
    int J2ab_t = state_Tz[ch_Tz][i].Jab_2;
    int a_t = state_Tz[ch_Tz][i].a;
    int b_t = state_Tz[ch_Tz][i].b;
    int c_t = state_Tz[ch_Tz][i].c;
    if (J2ab_t != J2ab || a_t != a)
      continue;
    if (b_t != b || c_t != c)
      continue;
    index = i;
    break;
  }
  if (index == -1) {
    cout << "wrong at find TP in TripleP::FindTP_Tz " << endl;
    cout << "J2 = " << J2 << "\t T_2z = " << T_2z << "\t ch_Tz = " << ch_Tz
         << "\t J2ab = " << J2ab << "\t a = " << a << "\t b = " << b
         << "\t c = " << c << endl;
    exit(0);
  }
  return index;
}

void TripleP::print() {
  for (int i = 0; i < channel.size(); i++) {
    cout << "===== JT channel : =====" << i << endl;
    cout << "T2 \t"
         << "T2ab \t"
         << "J2 \t"
         << "J2ab \t"
         << "a \t"
         << "b \t"
         << "c \t"
         << "E \t"
         << "index" << endl;
    for (int j = 0; j < state[i].size(); j++) {
      cout << state[i][j].T_2 << "\t";
      cout << state[i][j].Tab_2 << "\t";
      cout << state[i][j].J_2 << "\t";
      cout << state[i][j].Jab_2 << "\t";
      cout << state[i][j].a << "\t";
      cout << state[i][j].b << "\t";
      cout << state[i][j].c << "\t";
      cout << state[i][j].E << "\t";
      cout << state[i][j].index << "\t";
      cout << endl;
    }
  }
}

void TripleP::print(int channel_num) {
  // for(int i=0;i<channel.size();i++){
  int i = channel_num;
  cout << "===== JT channel : =====" << i << endl;
  cout << "T2 \t"
       << "T2ab \t"
       << "J2 \t"
       << "Jab \t"
       << "a \t"
       << "b \t"
       << "c \t"
       << "E \t"
       << "index" << endl;
  for (int j = 0; j < state[i].size(); j++) {
    cout << state[i][j].T_2 << "\t";
    cout << state[i][j].Tab_2 << "\t";
    cout << state[i][j].J_2 << "\t";
    cout << state[i][j].Jab_2 << "\t";
    cout << state[i][j].a << "\t";
    cout << state[i][j].b << "\t";
    cout << state[i][j].c << "\t";
    cout << state[i][j].E << "\t";
    cout << state[i][j].index << "\t";
    cout << endl;
  }
  //}
}

void TripleP::print(std::string filename) {
  std::ofstream file;
  file.open(filename);
  file << "a b c JJ-coupling state" << endl;
  file << "Energy cut-off is : " << Emax << " ( hbar-omega )" << endl;
  file << "total number of basis is : " << endl;
  int total_num = 0;
  for (int i = 0; i < state.size(); i++) {
    total_num += state[i].size();
  }
  file << total_num << endl;
  file << channel.size() << "\t channel number" << endl;
  for (int i = 0; i < channel.size(); i++) {
    file << state[i].size() << "\t";
  }
  file << endl;
  for (int i = 0; i < state.size(); i++) {
    file << "===== JT channel : =====" << i << endl;
    file << "T2 \t"
         << "T2ab \t"
         << "J2 \t"
         << "Jab \t"
         << "a \t"
         << "b \t"
         << "c \t"
         << "E \t"
         << "T2z \t"
         << "index" << endl;
    //<< "F \t"
    //<< "I \t" << endl;
    for (int j = 0; j < state[i].size(); j++) {
      file << state[i][j].T_2 << "\t";
      file << state[i][j].Tab_2 << "\t";
      file << state[i][j].J_2 << "\t";
      file << state[i][j].Jab_2 << "\t";
      file << state[i][j].a << "\t";
      file << state[i][j].b << "\t";
      file << state[i][j].c << "\t";
      file << state[i][j].E << "\t";
      file << state[i][j].T_2z << "\t";
      file << state[i][j].index << "\t";
      file << state[i][j].F << "\t";
      file << state[i][j].I << "\t";
      /*if (state[i][j].F == 1)
        file << "*";
      file << "\t";
      if (state[i][j].I == 1)
        file << "* \t";/**/
      file << endl;
    }
  }
  file.close();
}

void TripleP::print_Tz(std::string filename) {
  std::ofstream file;
  file.open(filename);
  file << "a b c JJ-coupling state_Tz" << endl;
  file << "Energy cut-off is : " << Emax << " ( hbar-omega )" << endl;
  file << "total number of basis is : " << endl;
  int total_num = 0;
  for (int i = 0; i < state_Tz.size(); i++) {
    total_num += state_Tz[i].size();
  }
  file << total_num << endl;
  file << channel_Tz.size() << "\t channel_Tz number" << endl;
  for (int i = 0; i < channel_Tz.size(); i++) {
    file << state_Tz[i].size() << "\t";
  }
  file << endl;
  for (int i = 0; i < channel_Tz.size(); i++) {
    file << "===== J channel_Tz : =====" << i << endl;
    file << "J2 \t"
         << "Jab \t"
         << "a \t"
         << "b \t"
         << "c \t"
         << "E \t"
         << "T2ab_z \t"
         << "T2z \t"
         << "index \t F\t I" << endl;
    for (int j = 0; j < state_Tz[i].size(); j++) {
      file << state_Tz[i][j].J_2 << "\t";
      file << state_Tz[i][j].Jab_2 << "\t";
      file << state_Tz[i][j].a << "\t";
      file << state_Tz[i][j].b << "\t";
      file << state_Tz[i][j].c << "\t";
      file << state_Tz[i][j].E << "\t";
      file << state_Tz[i][j].Tab_2z << "\t";
      file << state_Tz[i][j].T_2z << "\t";
      file << state_Tz[i][j].index << "\t";
      file << state_Tz[i][j].F << "\t";
      file << state_Tz[i][j].I << "\t";

      /*if (state_Tz[i][j].F > 0)
        file << "*";
      file << "\t";
      if (state_Tz[i][j].I > 0)
        file << "* \t";/**/

      file << endl;
    }
  }
  file.close();
}

void TripleP::print_TzM(std::string filename) {
  std::ofstream file;
  file.open(filename);
  file << "a b c m-scheme state_TzM" << endl;
  file << "Energy cut-off is : " << Emax << " ( hbar-omega )" << endl;
  file << "total number of basis is : " << endl;
  int total_num = 0;
  file << state_TzM.size() << endl;

  file << "T2z \t"
       << "J2m \t"
       << "a \t"
       << "b \t"
       << "c \t"
       << "T2ab_z \t"
       << "J2ab_m \t"
       << "E \t"
       << "index" << endl;
  for (int j = 0; j < state_TzM.size(); j++) {
    file << state_TzM[j].T_2z << "\t";
    file << state_TzM[j].J_2m << "\t";
    file << state_TzM[j].a_m << "\t";
    file << state_TzM[j].b_m << "\t";
    file << state_TzM[j].c_m << "\t";
    file << state_TzM[j].Tab_2z << "\t";
    file << state_TzM[j].Jab_2m << "\t";
    file << state_TzM[j].E << "\t";
    file << state_TzM[j].index << "\t";
    file << endl;
  }
  file.close();
}

void TripleP::print_channel() {
  int size = this->channel.size();
  cout << " TripleP channel Total Num: " << size << endl;
  cout << "index \t"
       << "num \t"
       << "J2 \t"
       << "T_2 \t"
       << "E_min \t"
       << "E_max \t" << endl;
  for (int i = 0; i < size; i++) {
    cout << i << "\t" << channel[i].num << "\t" << channel[i].J_2 << "\t"
         << channel[i].T_2 << "\t" << channel[i].E_min << "\t"
         << channel[i].E_max << endl;
  }
}
void TripleP::print_channel(std::string filename) {
  std::ofstream file;
  file.open(filename);
  int size = this->channel.size();
  file << " TP channel Total Num: " << size << endl;
  file << "index \t"
       << "num \t"
       << "J2 \t"
       << "T_2 \t"
       << "E_min \t"
       << "E_max \t" << endl;
  for (int i = 0; i < size; i++) {
    file << i << "\t" << channel[i].num << "\t" << channel[i].J_2 << "\t"
         << channel[i].T_2 << "\t" << channel[i].E_min << "\t"
         << channel[i].E_max << endl;
  }
  file << "E_min and E_max index:" << endl;
  file << "\t \t   "
       << "E : 0\t1\t2 \t3 \t4 \t5 ....." << endl;
  for (int i = 0; i < size; i++) {
    file << "channel : " << i << endl;
    file << "\t E_min index : ";
    for (int j = 0; j < channel[i].indexE_beg.size(); j++) {
      file << channel[i].indexE_beg[j] << "\t";
    }
    file << endl;
    file << "\t E_max index : ";
    for (int j = 0; j < channel[i].indexE.size(); j++) {
      file << channel[i].indexE[j] << "\t";
    }
    file << endl;
  }
  file.close();
}

void TripleP::print_channel_Tz(std::string filename) {
  std::ofstream file;
  file.open(filename);
  int size = this->channel_Tz.size();
  file << " TP channel_Tz Total Num: " << size << endl;
  file << "index \t"
       << "num \t"
       << "J2 \t"
       << "T_2z \t"
       << "E_min \t"
       << "E_max \t" << endl;
  for (int i = 0; i < size; i++) {
    file << i << "\t" << channel_Tz[i].num << "\t" << channel_Tz[i].J_2 << "\t"
         << channel_Tz[i].T_2z << "\t" << channel_Tz[i].E_min << "\t"
         << channel_Tz[i].E_max << endl;
  }
  file << "E_min and E_max index:" << endl;
  file << "\t \t   "
       << "E : 0\t1\t2 \t3 \t4 \t5 ....." << endl;
  for (int i = 0; i < size; i++) {
    file << "channel : " << i << endl;
    file << "\t E_min index : ";
    for (int j = 0; j < channel_Tz[i].indexE_beg.size(); j++) {
      file << channel_Tz[i].indexE_beg[j] << "\t";
    }
    file << endl;
    file << "\t E_max index : ";
    for (int j = 0; j < channel_Tz[i].indexE.size(); j++) {
      file << channel_Tz[i].indexE[j] << "\t";
    }
    file << endl;
  }
  file.close();
}
void TripleP::print_TPconf_info(std::string filename1, std::string filename2,
                                std::string filename3) {
  std::ofstream file;
  // file.open("data/TPconf.dat");
  file.open(filename1);
  int TPconf_size = 0;
  int TPconf_size_k = 0;
  int file_num = TPconf_k.size();
  for (int ch = 0; ch < this->TPconf.size(); ch++) {
    TPconf_size += this->TPconf[ch].size();
  }
  for (int k = 0; k < this->TPconf_k.size(); k++) {
    for (int ch = 0; ch < this->TPconf_k[k].size(); ch++) {
      TPconf_size_k += this->TPconf_k[k][ch].size();
    }
  }
  file << "TPconf_size  : " << TPconf_size
       << "\t TPconf_size_k : " << TPconf_size_k << endl;
  file << "File num : " << TPconf_k.size() << endl;
  for (int k = 0; k < TPconf_k.size(); k++) {
    file << "\t === file No. " << k << " ===" << endl;
    for (int ch = 0; ch < TPconf_k[k].size(); ch++) {
      file << "\t \t ch : " << ch << "\t size : " << TPconf_k[k][ch].size()
           << endl;
    }
  }
  for (int ch = 0; ch < TPconf.size(); ch++) {
    file << "\t == channel : " << ch << "\t == " << endl;
    for (int i = 0; i < TPconf[ch].size(); i++) {
      file << i << "\t" << TPconf[ch][i].F << "\t" << TPconf[ch][i].I << endl;
    }
  }
  file.close();
  // file.open("data/TPconf_E.dat");
  file.open(filename2);
  file << "TPconf_E : " << endl;
  for (int k = 0; k < file_num; k++) {
    for (int ch = 0; ch < TPconf_E[k].size(); ch++) {
      file << "\t ---- File_num : " << k << " ----- ";
      file << "\t == channel : " << ch << "\t == " << endl;
      for (int i = 0; i < TPconf_E[k][ch].size(); i++) {
        file << i << "\t" << TPconf_E[k][ch][i].E_F << "\t"
             << TPconf_E[k][ch][i].E_I << endl;
      }
    }
  }
  file.close();
  // file.open("data/TPconf_E_tot.dat");
  file.open(filename3);
  file << "TPconf_E_tot : " << endl;
  for (int k = 0; k < file_num; k++) {
    file << "\t ---- File_num : " << k << " ----- " << endl;
    for (int i = 0; i < TPconf_E_tot[k].size(); i++) {
      file << i << "\t" << TPconf_E_tot[k][i].E_F << "\t"
           << TPconf_E_tot[k][i].E_I << endl;
    }
  }
  file.close();
}

void TripleP::print_TPconf_info_Tz(std::string filename1,
                                   std::string filename2) {
  int TPconf_size_Tz = 0;
  int TPconf_size_k_Tz = 0;
  int file_num = TPconf_k_Tz.size();
  for (int ch = 0; ch < this->TPconf_Tz.size(); ch++) {
    TPconf_size_Tz += this->TPconf_Tz[ch].size();
  }
  for (int k = 0; k < this->TPconf_k_Tz.size(); k++) {
    for (int ch = 0; ch < this->TPconf_k_Tz[k].size(); ch++) {
      TPconf_size_k_Tz += this->TPconf_k_Tz[k][ch].size();
    }
  }

  std::ofstream file;
  // file.open("data/TPconf_Tz.dat");
  file.open(filename1);
  file << "TPconf_size_Tz : " << TPconf_size_Tz
       << "\t TPconf_size_k_Tz : " << TPconf_size_k_Tz << endl;
  for (int k = 0; k < TPconf_k_Tz.size(); k++) {
    file << "\t === file No. " << k << " ===" << endl;
    for (int ch = 0; ch < TPconf_k_Tz[k].size(); ch++) {
      file << "\t \t ch : " << ch << "\t size : " << TPconf_k_Tz[k][ch].size()
           << endl;
    }
  }

  for (int ch = 0; ch < TPconf_Tz.size(); ch++) {
    file << "\t == channel : " << ch << "\t == " << endl;
    for (int i = 0; i < TPconf_Tz[ch].size(); i++) {
      file << i << "\t" << TPconf_Tz[ch][i].F << "\t" << TPconf_Tz[ch][i].I
           << endl;
    }
  }
  file.close();

  // file.open("data/TPconf_E_Tz.dat");
  file.open(filename1);
  file << "TPconf_E_Tz : " << endl;
  for (int k = 0; k < file_num; k++) {
    for (int ch = 0; ch < TPconf_E_Tz[k].size(); ch++) {
      file << "\t ---- File_num : " << k << " ----- ";
      file << "\t == channel : " << ch << "\t == " << endl;
      for (int i = 0; i < TPconf_E_Tz[k][ch].size(); i++) {
        file << i << "\t" << TPconf_E_Tz[k][ch][i].E_F << "\t"
             << TPconf_E_Tz[k][ch][i].E_I << endl;
      }
    }
  }
  file.close();
}

int TripleP::Tri(int &j2a, int &j2b, int &j2c) {
  // I+J>=K, I+K>=J, J+K>=I,
  // I/2+J/2+K/2 = INTEGER.
  // TRI=1, WHEN TRIADIC CONDITION IS FULFILLED, TRI=-1 OTHERWISE
  int L2 = j2a + j2b + j2c;
  if (int(L2) != (L2 / 2) * 2)
    return -1;
  L2 = L2 / 2;
  if (j2a * j2b * j2c < 0)
    return -1;
  if ((L2 - j2a) * (L2 - j2b) * (L2 - j2c) < 0)
    return -1;
  return 1;
}

int TripleP::TzDiff(State_TP &A, State_TP &B) {
  int res = -1;
  if (A.a != B.a || A.b != B.b || A.c != B.c)
    return res;
  if (A.Jab_2 != B.Jab_2 || A.Tab_2 != B.Tab_2)
    return res;
  if (A.J_2 != B.J_2 || A.T_2 != B.T_2 || A.E != B.E)
    return res;
  res = 1;
  return res;
}
//~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~end Triple
// Particle~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~//

//~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~end
// Check~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~//
